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ABSTRACT 

This  paper  summarizes  a  new  methodology  to  design  sequential  non-diagonal  QFT  controllers  for  multi- 
input-multi-output  MIMO  systems  with  uncertainty,  which  is  a  central  issue  in  UAV  control  systems.  It 
also  demonstrates  the  feasibility  of  that  methodology  to  control  the  position  and  attitude  of  a  6x6  MIMO 
spacecraft  with  large  flexible  appendages.  The  last  part  of  the  paper  introduces  a  new  practical 
methodology  to  design  robust  controllers  that  work  under  a  switching  mechanism,  going  beyond  the 
classical  linear  limitations  and  giving  a  solution  for  the  well-known  robustness-performance  trade-off. 

1.0  INTRODUCTION 

Control  of  multivariable  systems  (multiple-input-multiple-output,  MIMO)  with  model  uncertainty  is  still 
one  of  the  hardest  problems  that  control  engineers  have  to  face  in  Unmanned  Air  Vehicle  (UAV)  real- 
world  applications.  Input-output  directionality,  coupling  among  control  loops,  transmission  zeros,  pairing, 
etc.  are  some  of  the  main  complexities  that  define  a  MIMO  system.  Moreover,  model  uncertainties 
substantially  increase  such  difficulties,  making  more  restrictive  the  inherent  performance  limitations  of  the 
control  system.  In  the  last  few  decades  a  very  significant  amount  of  work  in  linear  MIMO  systems  has 
been  done.  The  first  technique  that  made  a  quantitative  synthesis  of  MIMO  systems,  taking  into  account 
quantitative  bounds  on  the  plant  uncertainty  and  quantitative  tolerances  on  the  acceptable  closed-loop 
system  response,  was  the  Quantitative  Feedback  Theory  (QFT)  [1].  In  the  last  few  years  some  new 
methods  for  non-diagonal  (full  matrix)  multivariable  QFT  robust  control  system  design  have  been 
introduced.  The  first  part  of  the  paper  introduces  a  new  methodology  [2-6]  that  improves  the  current  non¬ 
diagonal  MIMO  QFT  control  techniques.  The  second  part  validates  the  new  techniques  by  applying  them 
to  control  the  position  and  attitude  of  a  6x6  spacecraft  with  large  flimsy  appendages  [7]. 

Combining  robust  designs  and  stable  switching,  the  control  strategy  could  optimize  the  time  response  of 
the  system  by  fast  adaptation  of  the  controller  parameters  during  the  transient  response  according  to 
certain  rules  based  on  the  amplitude  of  the  error.  The  last  part  of  the  paper  introduces  a  methodology  to 
design  a  family  of  robust  controllers  able  to  go  beyond  the  classical  linear  performance  limitations.  The 
methodology  is  based  on  both  a  new  graphical  stability  criterion  for  switching  linear  systems  and  the 
robust  quantitative  feedback  theory  (QFT)  [8]. 

2.0  NON-DIAGONAL  MIMO  QFT  CONTROL  DESIGN  METHODOLOGY  [2  7] 

Control  of  multivariable  systems  (multiple-input-multiple-output,  MIMO)  with  model  uncertainty  are  still 
one  of  the  hardest  problems  that  the  control  engineer  has  to  face  in  real-world  applications.  Three  of  the 
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main  characteristics  that  define  a  MIMO  system  are  the  input  and  output  directionality  -different  vectors  to 
actuate  U  and  to  measure  F-;  the  coupling  among  control  loops  -some  outputs  y-x  can  be  influenced  by 
several  inputs  wi?  and  some  inputs  ux  can  influence  several  outputs  y{;  and  the  transmission  zeros  of  the 
plant  matrix. 

In  the  last  few  decades  a  very  significant  amount  of  work  in  MIMO  systems,  too  numerous  to  list  here,  has 
been  done.  Using  MIMO  QFT,  Horowitz  proposed  to  translate  the  original  nxn  MIMO  problem  into  n 
separate  quantitative  multiple-input-single-output  MISO  problems,  each  with  plant  uncertainty,  external 
disturbances  and  closed-loop  tolerances  derived  from  the  original  problem  [1].  Two  different  approaches, 
the  so-called  sequential  and  non-sequential  methods,  consider  in  successive  iterative  steps  an  equivalent 
plant  that  either  takes  also  into  account  the  controllers  designed  in  the  previous  steps,  or  only  deals  with 
the  plant  respectively. 

However,  although  such  original  MIMO  QFT  methods  take  the  coupling  among  loops  into  account,  they 
only  propose  the  use  of  a  diagonal  controller  G  to  govern  the  MIMO  plant.  This  structure  can  be  improved 
using  non-diagonal  controllers.  In  fact,  a  fully  populated  matrix  controller  allows  the  designer  much  more 
design  flexibility  to  control  MIMO  plants  than  the  classical  diagonal  controller  structure.  The  use  of  the 
non-diagonal  components  can  also  ease  the  diagonal  controller  design  problem.  In  the  last  few  years  some 
new  methods  for  non-diagonal  multivariable  QFT  robust  control  system  design  have  been  introduced.  For 
the  sake  of  clarity,  this  section  summarizes  a  previous  work  [2-7]  that  extends  the  classical  QFT  diagonal 
controller  design  for  MIMO  plants  with  uncertainty  to  the  fully  populated  matrix  controller  design.  The 
work  studies  three  cases:  the  reference  tracking,  the  external  disturbance  rejection  at  plant  input  and  the 
external  disturbance  rejection  at  plant  output.  It  presents  the  definition  of  three  specific  coupling  matrices 
(ciij5  c2ij,  C3ij),  one  for  each  case,  and  introduces  a  sequential  design  methodology  for  non-diagonal  QFT 
controllers. 


2.1  The  Coupling  Matrix 

The  objective  of  this  section  is  to  define  a  measurement  index  (the  coupling  matrix)  that  allows  one  to 
quantify  the  loop  interaction  in  MIMO  control  systems.  Consider  a  nxn  linear  multivariable  system  -see 
Fig.  1-,  composed  of  a  plant  P,  a  fully  populated  matrix  controller  G,  a  pre-filter  F,  a  plant  input 
disturbance  transfer  function  Pd i?  and  a  plant  output  disturbance  transfer  function  Pdo,  where  P  e  3P  ,  3P 
is  the  set  of  possible  plants  due  to  uncertainty,  and, 


H*)= 

Pi  i(s)  Pn(s)  ...  Puis) 
P2\(S)  P2l{s)  -  Pin  ('*’) 

;  G(s)  = 

g  n(s)  gi2(s)  -  gjs) 
g2l(s)  gn(s)  -  &2„(s') 

;  F(,)  = 

fn{s)  fn{s) 
fl\{S)  fn  CO  •• 

••  AM 

-  fln{s) 

_Pn\iS)  Pnl(s)  -  Pnn{S)_ 

_g„lG')  gM  -  gm{s)_ 

_fn\(S)  fnl{S)  •' 

-  fnn(S)_ 

(1) 

The  reference  vector  r  ’  and  the  external  disturbance  vectors  at  plant  input  '  and  plant  output  d0  ’  are  the 
inputs  of  the  system.  The  output  vector  y  is  the  variable  to  be  controlled. 

It  is  denoted  P  as  the  plant  inverse  so  that, 


p{sy=P*(s)=[P:{s)\=A(s)+B{s) 


- 1 

S3* 

0 

0 

0 

...  /^(s) 

0 

0 

+ 

0 

0 

0 

P*nn(S)  _ 

/nAS) 

0 

(2) 
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G(,)=Gd(.S)+Gb(.s) 


Sn(s) 

o 

o 

j 

0 

•••  gl  n(S) 

0 

0 

+ 

0 

0 

0  Snn(s)_ 

_g«l(s) 

0 

(3) 


where  A  is  the  diagonal  part  and  B  is  the  balance  of  P*;  and  Gd  is  the  diagonal  part  and  Gb  is  the  balance  of 
G.  The  next  paragraphs  introduce  a  measurement  index  to  quantify  the  loop  interaction  in  the  three 
classical  cases:  reference  tracking,  external  disturbances  at  plant  input,  and  external  disturbances  at  plant 
output.  That  index  is  called  the  coupling  matrix  and,  depending  on  the  case,  shows  three  different 
expressions:  C\,  C2,  C3  respectively. 


Fig.  1  Structure  of  a  2  Degree  of  Freedom  MIMO  System 


2.1.1  Tracking 

The  transfer  function  matrix  of  the  controlled  system  for  the  reference  tracking  problem,  without  any 
external  disturbance,  can  be  written  as  shown  in  Eq.  (4), 


y  =  (I  +  PG)  1  P  G  r  =  Ty  /•  =  Ty  F  /•'  (4) 

Using  Eq.  (2)  and  (3),  Eq.  (4)  can  be  rewritten  as, 

Ty/r  r  =  {l  +  A'1  (jd)-1  A'1  Gd  r  +  il  +  A-1  Gd)~l  A'1  {Gb  r-(B  +  Gh)  Ty/r  r )  (5) 

In  the  expression  of  the  closed-loop  transfer  function  matrix  of  Eq.  (5),  it  is  possible  to  find  two  different 
terms: 

i.  A  diagonal  term  Ty/r_ d, 

Ty/r_d  =  {I  +  G& )  Gd  (6) 

that  presents  a  diagonal  structure.  Note  that  it  does  not  depend  on  the  non-diagonal  part  of  the  plant 
inverse  B,  nor  on  the  non-diagonal  part  of  the  controller  It  is  equivalent  to  n  reference  tracking  SISO 
systems  formed  by  plants  equal  to  the  elements  of  A'1  when  the  n  corresponding  parts  of  a  diagonal  Gd 
control  them,  as  shown  in  Fig.  2a. 
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ii.  A  non-diagonal  term  Jy/r_ b, 

Ty,r_ b  =  (/  +  A-1  Gd  )“'  [<?b -(B  +  Gb)  Ty/r  ]=(/  +  A-1  Gd  f  /l"1  C,  (7) 

that  presents  a  non-diagonal  structure.  It  is  equivalent  to  the  same  n  previous  systems  with  internal 
disturbances  cx-  at  plant  input  (Fig.  2b). 

In  Eq.  (7),  the  matrix  C\  is  the  only  part  that  depends  on  the  non-diagonal  parts  of  both  the  plant  inverse  B 
and  the  controller  G  b  .  Hence,  it  comprises  the  coupling,  and  from  now  on  C\  will  be  the  coupling  matrix 
of  the  equivalent  system  for  reference  tracking  problems, 

Ci  =  Gb  -  (B  +  Gb )  Ty/r  (8) 

Each  element  cVx]  of  this  matrix  obeys, 

n 

clij  =  Sij  0  -  §ij )-  X  +  £ik  )  0  -  5ik )  (9) 

k  =  1 

where  5ki  is  the  delta  of  Kronecker  that  is  defined  as, 


|4  =  1<=>k  =  i 

[4,  =0<»k^i 


(a) 

n 


(10) 


Fig.  2  i-th  equivalent  SISO  and  MISO  systems 


2.1.2  Disturbance  rejection  at  plant  input 

The  transfer  matrix  from  the  external  disturbance  at  plant  input  d\  to  the  output  y  can  be  written  as  shown 
in  Eq.  (11), 
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y  =  (l  +  P  G)-]  P  di  =  Ty/di  dt  =  Tyldi  Pdi  d\  (11) 

and  then, 

T„  d,  =[l  +  A  '  a,)'  A-'  d  -(l  +  A'  a,}'  A'  (( B+G„)TyJd ,  (12) 

In  that  expression  -Eq.  (12)-  it  is  possible  to  find  two  different  terms: 

i.  A  diagonal  term  Ty/ di_d, 

Tyj«j=(l  +  A'1  <?„('  A'1  (13) 

Again,  Eq.  (13)  is  equivalent  to  n  regulator  MISO  systems,  as  shown  in  Fig.  3a. 

ii.  Non  diagonal  term  Jy/di_b 

V »  =  (/  +  A''  G„ A-'  (B+  C„ ) =(l  +  X'  GX  A'  C2  (14) 

that  presents  a  non-diagonal  structure  which  is  equivalent  to  the  same  n  previous  systems  with  external 
disturbances  c2ij-  di}  at  plant  input,  as  shown  in  Fig.  3b. 


In  Eq.  (14),  the  matrix  C2  comprises  the  coupling,  and  from  now  on  C2  will  be  the  coupling  matrix  of  the 
equivalent  system  for  external  disturbance  rejection  at  plant  input  problems, 

C2=(B  +  Gb)Ty/di  (15) 


(a) 

n 


Fig.  3  i-th  equivalent  MISO  systems 
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Each  element  c2\}  of  this  matrix  obeys, 

c2ij  =  Z  ( Ak  +  Ak )  4j  0  -  4 )  (16> 

k=l 

where  5ki  is  the  delta  of  Kronecker  defined  in  Equation  (10). 

2.1.3  Disturbance  rejection  at  plant  output 

The  transfer  matrix  from  the  external  disturbance  at  plant  output  d0  to  the  output  y  can  be  written  as 
shown  in  Eq.  (17), 

y  =  (I  +  P  G)~l  d0  =Ty/dod0=  Ty/do  Pdo  d0  (17) 

and  then, 

V»rf»=(/  +  /l'lcj)"1  da+(l  +  A-'  Gdy'  A-'  {B-(B+Gb)Ty/Jd0  (18) 

In  that  expression  -Eq.  (18)-  it  is  possible  to  find  two  different  terms: 

i.  A  diagonal  term  Jy/d0_d? 

Ty/do_d  =(l  +  Gd  )  (19) 

Once  more,  Eq.  (19)  is  equivalent  to  the  n  regulator  MISO  systems  showed  in  Fig.  4a, 

ii.  Non  diagonal  term  Jy/do_b 

ry/d0_b  =  (/  +  A-1  Gd  Y  A-1  [b  -  {B  +  Gb )  Ty/io  ]  =  (/  +  A-1  Gd  f  A1  C3  (20) 

that  presents  a  non-diagonal  structure.  It  is  equivalent  to  the  same  n  previous  systems  with  external 
disturbances  c3i j  doj  at  plant  input,  as  shown  Fig.  4b. 

In  Eq.  (20),  the  matrix  C3  comprises  the  coupling,  and  from  now  on  it  will  be  the  coupling  matrix  of  the 
equivalent  system  for  external  disturbance  rejection  at  plant  output  problems, 

C3=B-(B  +  Gh )  Ty/do  (21) 

Each  element  of  the  coupling  matrix,  c3ij  obeys, 

c3ij  =Pn  (l-^j)-Z ( Ak  +  gik )  4j  (!  - 4 )  (22) 

k= 1 

where  5ki  is  the  delta  of  Kronecker  as  defined  in  Equation  (10). 
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(a) 


n 


Fig.  4  i-th  equivalent  MISO  systems 


2.2  The  Coupling  Elements 

In  order  to  design  a  MIMO  controller  with  a  low  coupling  level,  it  is  necessary  to  study  the  influence  of 
every  non-diagonal  element  gy  on  the  coupling  elements  c\\ j,  c2y  and  c3ij,  defined  in  Eq.  (9),  (16)  and  (22). 
These  elements  can  be  simplified  to  quantify  the  coupling  effects.  Then  it  will  be  possible  to  analyze  the 
loop  decoupling  and  to  state  some  conditions  and  limitations  using  fully  populated  matrix  controllers.  To 
analyze  the  coupling  elements,  one  Hypothesis  is  stated. 


Hypothesis  HI:  suppose  that  in  Eq.  (9),  (16)  and  (22), 


/  *  \ 

/  *  \ 

y^ij  +  g'ijf  jj 

» 

Wik"l"  <?ik  kj 

for  k  ^  j,  and  in  the  bandwidth  of  t  jj 


(23) 


Note  that  the  above  expression  is  scale  invariant  and  is  typically  fulfilled  once  the  MIMO  system  has  been 
ordered  according  to  appropriate  methods  like  the  Relative  Gain  Analysis,  etc.  Then  the  diagonal  elements 
/jj  will  be  much  larger  that  the  non-diagonal  ones  /kj, 

tjj|»  U  for  k  ^  j,  and  in  the  bandwidth  of  (24) 

Now,  two  simplifications  are  applied  to  facilitate  the  quantification  of  the  coupling  effects  chj,  c2 ij,  c3ij. 

Simplification  SI:  Using  the  Hypothesis  HI,  Eqs.  (9),  (16)  and  (22),  which  describe  the  coupling  elements 
in  the  tracking  problem,  disturbance  rejection  at  plant  input  and  disturbance  rejection  at  plant  output 
respectively,  are  rewritten  as  shown  Table  I. 

Simplification  S2:  The  elements  t jj  are  computed  for  each  case  from  the  equivalent  system  derived  from 
Eqs.  (6),  (13)  and  (19).  The  results  are  shown  in  Table  I. 
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Table  I.  Simplifications  to  quantify  the  coupling  effects 


Reference  tracking 

External  disturbances 
at  plant  input 

External  disturbances 
at  plant  output 

Simplification 

SI 

ciij=£ij-'jj  Ui+gij) ;  i^j 

(25) 

c2ij=tjj  Ui+gij)  ;  i^J 

(26) 

C3ij=Aj  -'v  Uj+gi;)  ;  j 

(27) 

Simplification 

S2 

<i=  s" m 

i  +  gjj  Pjj 

<»  =  p"  <*>> 
l  +  Sa  Pa 

II 

+ 

Op 

IS 

o 

Due  to  Simplifications  SI  and  S2,  the  coupling  effects  c hj,  c2 ij,  c3 y  can  be  computed  as, 
Tracking 


clij=£ij- 


gjj  faj+gjj) 

(pii+gji) 


i*j 


Disturbance  rejection  at  plant  input 

fei+^ij)  . 


c2ij 


(Pjj  +  gjj) 


Disturbance  rejection  at  plant  output 


C3ij 


* 

Pa- 


*  /  *  \ 

Pti  bii  +  gij) 

IpI+Su) 


i^j 


(31) 


(32) 


(33) 


2.3  The  Optimum  Non-diagonal  Controller 

Consider  non-diagonal  controllers  to  reduce  the  coupling  effect  and  diagonal  controllers  that  help  to 
achieve  the  loop  performance  specifications.  The  optimum  non-diagonal  controllers  for  the  three  cases 
(tracking  and  disturbance  rejection  at  plant  input  and  output)  can  be  obtained  making  the  loop  interaction 
of  Eqs.  (31),  (32)  and  (33)  equal  to  zero. 

Note  that  both  elements,  p*-}  and  • ,  of  these  equations  are  uncertain  elements  of  P*.  Every  uncertain 
plant  /?ij  can  be  any  plant  represented  by  the  family, 

{Pij}=.PijN(1  +  Aij)  >  0 -|Au|  -  *  fori,j  =  l,...,n  (34) 


where 


Pii 


is  the  nominal  plant,  and  A/?*-  the  maximum  of  the  non-parametric  uncertainty  radii 


The  nominal  plants 


*N 

Pii 


J  *N 
and  p  jj 


that  will  be  chosen  for  the  optimum  non-diagonal  controller  will 
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follow  the  next  rules: 

a)  If  the  uncertain  parameters  of  the  plants  show  a  uniform  Probability  Distribution  (Fig.  5a)  -which  is 

*  * 

typical  in  the  QFT  methodology-,  then  the  elements  and  p jj  for  the  optimum  non-diagonal 

controller  will  be  the  nominal  plants  p y  and  p^  ,  which  minimise  the  maximum  of  the  non- 

parametric  uncertainty  radii  A p*-  and  A p*-  that  comprise  the  plant  templates  (Fig.  5b). 

b)  If  the  uncertain  parameters  of  the  plants  show  a  non-uniform  Probability  Distribution  (Fig.  5c),  then 

*  *  *N 

the  elements  p{j  and  p jj  for  the  optimum  non-diagonal  controller  will  be  the  nominal  plants  px j 

*N 

and  p  jj  ,  whose  set  of  parameters  maximize  the  area  of  the  Probability  Distribution  in  the  regions 

[  a{-  -  £ ,  a{-  +  £  ]  and  [a---£,a--+s](\/  parameter  a^9  . . .,  tty,  . . .)  respectively. 

Now,  making  Eqs.  (31),  (32)  and  (33)  equal  to  zero  and  using  Eq.  (34),  the  optimum  non-diagonal 
controller  for  each  case  is  obtained. 


Uniform 

PD 


(a) 


Parameter  ay 


! l 

aij_min  ^ij_max 


Fig.  5  Probability  Distribution  of  the  parameter  ay,  and  Non-parametric  uncertainty  radii  Ap^  that  comprise 

the  plant  templates 
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2.3.1  Tracking 


sT=Fpd 


*ji 


* 


*N 

Pjj  , 


,  for  i  ^  j 


2.3.2  Disturbance  rejection  at  plant  input 


sT=Fpd 


* 


,  for  i  ^  j 


2.3.3  Disturbance  rejection  at  plant  output 


gT=Fp“ 


£jj 


*N 
*  N 


,  for  i  ^  j 


(35) 


(36) 


(37) 


l  P  jJ  J 

where  the  function  Fpd(A )  means  in  every  case  a  casual  and  stable  proper  function  made  from  the 
dominant  poles  and  zeros  of  the  expression  A. 


2.4  The  Coupling  Effects 

The  minimum  achievable  coupling  effects  -Eqs.  (38),  (40),  (42)-  can  be  computed  substituting  the 
optimum  controller  of  Eqs.  (35),  (36)  and  (37)  in  the  coupling  expressions  of  Eqs.  (31),  (32)  and  (33) 
respectively,  and  taking  into  account  the  uncertainty  radii  of  Eq.  (34).  Analogously,  the  maximum 
coupling  effect  without  any  non-diagonal  controller  -pure  diagonal  controller  cases-  can  be  computed 
substituting  gij=0  in  the  Eqs.  (31),  (32)  and  (33)  respectively  -Eqs.  (39),  (41),  (43)-.  That  is  to  say, 


2.4.1 


Tracking 


=  ^i.i  (Aij  -  Ayh'y 

hi  (1  +  Aij)^jj| 


2.4.2 


Disturbance  rejection  at  plant  input 


2.4.3  Disturbance  rejection  at  plant  output 

=hK-Au)*j 

=  k  (1+Ahd 


r3ijL..=g°pt 

5  li  5i, 


Cv- 

3ij 


gij=° 


(38) 

(39) 


(40) 

(41) 


(42) 

(43) 


RTO-EN-SCI-195 


1  - 10 


NATO 

OTAN 


Beyond  the  Classical  Performance  Limitations  Controlling  Uncertain  MIMO  Systems 


ORGANIZATION 


where, 


and  the  uncertainty  is:  0  <  A-  |  <  Ap*-  ,  0  <  <  A/?*.  ,  for  i,  j  =  l,...,n 


(44) 


The  coupling  effects,  calculated  in  the  pure  diagonal  controller  cases,  result  in  three  expressions  (39),  (41) 
and  (43)  that  still  present  a  non-zero  value  when  the  nominal-actual  plant  mismatching  due  to  the 
uncertainty  disappears:  A-  =  0  and  A  -  =  0  .  However,  the  coupling  effects  obtained  with  the  optimum  non¬ 
diagonal  controllers  -Eqs.  (38),  (40)  and  (42)-  tends  to  zero  when  that  mismatching  disappears. 


2.5  Design  Methodology 

The  proposed  controller  design  methodology  is  a  sequential  procedure  closing  loops  with  four  steps  [2-7]: 

Step  A:  Controller  structure,  input-output  pairing  and  loop  ordering.  First,  the  methodology  identifies 
the  controller  structure  (minimum  required  elements  of  the  controller  matrix)  and  the  input-output  pairings 
by  using  the  frequency-dependent  Relative  Gain  Array  -RGA-  [10-11].  Then,  the  matrix  P  (s)  is 
reorganized  so  that  \pn  (s)]_1  has  the  smallest  phase  margin  frequency,  \p2 2  (s)]-1  the  next  smallest  phase 
margin  frequency,  and  so  on  to  guarantee  the  existence  of  a  solution  [1]. 


After  that,  the  sequential  design  technique  composed  of  n  stages,  as  many  as  loops,  performs  the  following 
two  steps  B  and  C  for  every  column  of  the  matrix  compensator  G(s )  from  k  =  1  to  n  (Fig.  6). 
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n 

Fig.  6  Steps  for  controllers  design 


Step  B:  Design  of  the  diagonal  compensator  gufs).  The  diagonal  element  g\±(s)  is  calculated  through 
standard  QFT  loop-shaping  [1]  for  the  inverse  of  the  equivalent  plant  [pkk  e  (s)]kl  in  order  to  achieve 
robust  stability  and  robust  performance  specifications  [13-14].  The  equivalent  plant  satisfies  the  recursive 
relationship  (45)  [13],  which  is  an  extension  for  the  non-diagonal  case  of  the  recursive  expression 
proposed  by  Horowitz  [12]  as  the  Improved  design  technique ,  also  called  Second  method  by  Houpis  et  al. 
[1]. 


r  .  1  r  ..  i  ( l/’ in:. 1, (A  ,  +  U,ik.i;i(''l],_i  HIamhU'IJ,.  !  +  U'-k-i.ji.v)],.-,) 

w 'ML  -  I*  ML - ,  .  J  , - -j - • 


(45) 


i, j  —  k;  [R*e(s)]k=1=R*(s) 
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If  the  control  system  requires  tracking  specifications  as 

t^x  “^rii  +  ^ciii"Eq*(5)">  the  tracking  bounds  bn  and  an  will 
specification  rc m,  so  that: 


aii(®)  -  (jco)  <  ^ii(co)  then,  because 
have  to  be  corrected  with  the  coupling 


b[[  b[[  -  Tc\\[  ,  Cl[[  ^ii^~^clii 

*clii  =  Wii  Clii  ^  Tcln 


(46) 

(47) 

(48) 


These  are  the  same  corrections  proposed  originally  by  Horowitz  (see  also  [1]).  However,  with  the 
proposed  non-diagonal  method  these  corrections  will  be  less  demanding.  The  coupling  expression  tcm  =  wn 
cm  is  now  minor  than  in  the  previous  diagonal  methods  -compare  Eqs.  (38)  and  (39)-.  The  off-diagonal 
elements  gy  (i^j)  of  the  matrix  controller  will  attenuate  or  cancel  that  cross  coupling.  Then  the  diagonal 
elements  of  the  non-diagonal  method  will  need  less  bandwidth  than  the  diagonal  elements  of  the 
previous  diagonal  methods. 

Step  C:  Design  of  the  (n-1)  non-diagonal  elements  gifs)  (i  ^ k ,  i  =  l,2,...n).  The  gifs)  (i  ^  k)  elements  of 
the  k-th  compensator  column  are  designed  to  minimize  the  non-diagonal  elements  of  the  cross-coupling 
matrices  according  to  different  purposes:  reference  tracking  (31),  (35);  disturbance  rejection  at  plant  input 
(32),  (36);  and  disturbance  rejection  at  plant  output  (33),  (37).  The  resulting  compensators  gk(s)  have  to  be 
casual  and  stable,  and  include  the  dominant  dynamics. 

The  off-diagonal  controller  elements  can  be  allocated  not  only  to  reduce  the  coupling  effects  of  the  MIMO 
system,  but  also  to  reach  complementary  objectives,  such  as  to  remove  RHP  (right-half  plane) 
transmission  zeros  introduced  during  the  controller  design  [5],  improve  system  integrity  [13]  and  stability 
margins,  reduce  controller  efforts,  etc. 

Step  D:  Design  of  the  prefilter.  The  design  of  the  prefilter  F(s)  does  not  present  any  additional  difficulty 
because  the  final  transfer  function  that  relates  R(s)  to  F(s)  shows  less  loop  interaction  thanks  to  the  fully 
populated  compensator  design.  Therefore,  the  prefilter  F(s)  can  generally  be  a  diagonal  matrix. 


2.6  Stability  Conditions 

Closed-loop  stability  of  a  MIMO  system  with  a  non-diagonal  controller  designed  by  using  a  sequential 
procedure  is  guaranteed  by  the  following  sufficient  conditions  [14]: 

(c.l)  each  L[(s)  =  gfs)  [pii*e(s)]f  \  i=l,  ...,  n ,  satisfies  the  Nyquist  encirclement  condition, 

(c.2)  no  RHP  pole-zero  cancellations  occur  between  gfs)  and  [/?;;*%?)],  ',  i=l 
(c.3)  no  Smith-McMillan  pole-zero  cancellations  occur  between  P(s)  and  G(s ),  and 
(c.4)  no  Smith-McMillan  pole-zero  cancellations  occur  in  |  P*(s)  +  G(s)  \ . 

2.7  Remarks 

It  is  important  to  note  that  the  calculation  of  the  equivalent  plant  [pkk*e(^)]k_1?  (45),  usually  introduces  some 
exact  pole-zero  cancellations.  That  operation  could  be  precisely  performed  by  using  symbolic 
mathematical  tools  [1].  However,  fictitious  poles  and  zeros  may  be  introduced  when  using  numerical 
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calculus  due  to  the  typical  rounding  errors  of  the  computer.  Additionally,  it  is  needed  to  determine  the 
inverse  of  the  plant  matrix,  which  can  also  be  numerically  non-reliable. 

In  this  paper,  these  problems  are  overcome  through  a  new  frequency  response  computation  method.  That 
is,  for  each  frequency  of  interest  co  and  for  every  set  of  parameters  within  the  region  of  uncertainty,  each 
element  Pi](jco)  of  the  plant  transfer  function  matrix  is  translated  into  a  complex  matrix  Pfreqj j  that 
represents  the  frequency  response  of  every  plant  element  within  the  uncertainty.  Thus,  this  complex  matrix 
has  as  many  rows  as  different  cases  generated  due  to  the  uncertainty  and  as  many  columns  as  frequencies 
(49).  All  the  abovementioned  calculations  are  then  performed  on  the  basis  of  this  set  of  complex  matrices 
by  using  element-by-element  matrix  operations.  As  a  result,  potential  impediments  related  to  practical 
computation  are  avoided. 


frequencies  =^> 

cox  co 2 

cor 

••• 

Case  1 

CL  i  x  $12 

•••  a\n 

Case  2 

>  _ 

a2\ 

/ra?Jj  Case  k 

Case  m 

_am\  . 1 

•  •  •  a 

mn 

"►  Plant  k 


akr  =  Rc(afo. )  +  j  Tmag(«/cr ) 


Template  at 


(49) 


At  the  same  time,  arbitrarily  picking  the  wrong  order  of  the  loops  to  be  designed  can  result  in  the  non¬ 
existence  of  a  solution.  This  may  occur  if  the  solution  process  is  based  on  satisfying  an  upper  limit  of  the 
phase  margin  frequency  co^  for  each  loop.  Hence,  Loop  i  having  the  smallest  phase  margin  frequency  will 
have  to  be  chosen  as  the  first  loop  to  be  designed.  The  loop  that  has  the  next  smallest  phase  margin 
frequency  will  be  next,  and  so  on  [1]. 

Although  very  remote,  theoretically  there  exists  the  possibility  of  introducing  RHP  transmission  zeros  due 
to  the  compensator  design.  This  undesirable  situation  can  not  be  detected  until  the  multivariable  system 
design  is  completed.  To  avoid  it  the  proposed  methodology  (Steps  A ,  B  and  Q  is  inserted  in  a  procedure 
introduced  by  Garcia-Sanz  and  Eguinoa  [5].  Once  the  matrix  compensator  G(s )  is  designed,  the 
transmission  zeros  of  P(s )  G(s )  are  determined  using  the  Smith-McMillan  form  and  over  the  set  of 
possible  plants  3P  due  to  uncertainty.  If  there  exist  new  RHP  transmission  zeros  apart  from  those  initially 
present  in  P(s ),  they  can  be  removed  by  using  the  non-diagonal  elements  placed  in  the  last  column  of  the 
matrix  G(s). 


3.0  MIMO  QFT  CONTROL  FOR  A  SPACECRAFT  WITH  LARGE  FLEXIBLE 
APPENDAGES  [7] 

This  section  summarizes  the  design  of  a  robust  non-diagonal  MIMO  QFT  controller  to  govern  the  position 
and  attitude  of  a  Darwin-type  spacecraft  with  large  flexible  appendages.  The  satellite  is  one  of  the  flyers  of 
a  multiple  spacecraft  constellation  for  a  future  ESA  mission.  It  presents  a  6x6  high  order  MIMO  model 
with  large  uncertainty  and  loop  interactions  introduced  by  the  flexible  modes  of  the  low-stiffness 
appendages.  The  scientific  objectives  of  the  satellite  require  very  demanding  control  specifications  for 
position  and  attitude  accuracy,  high  disturbance  rejection,  loop-coupling  attenuation  and  low  order 
controller.  This  section  demonstrates  the  feasibility  of  sequential  non-diagonal  MIMO  QFT  strategies 
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controlling  the  Darwin  spacecraft  and  compares  the  results  with  a  previous  H-infinity  design. 

3.1  Description 

The  Darwin  mission  consists  of  three  to  six  telescopes  arranged  in  a  symmetric  configuration  flying  in 
formation  around  a  master  satellite  or  central  hub  (Fig.  7).  Darwin  will  employ  nulling  interferometry  to 
detect  and  analyze  through  appropriate  spectroscopy  techniques  the  atmosphere  of  remote  planets  close  to 
a  bright  star.  The  infrared  light  collected  by  the  free  flying  telescopes  will  be  recombined  inside  the  hub- 
satellite  in  such  a  way  that  the  light  from  the  central  star  suffers  destructive  interference  and  is  cancelled 
out,  allowing  this  way  the  much  fainter  planet  easier  to  stand  out.  The  interferometry  requires  very 
accurate  and  stable  positioning  of  the  spacecraft  in  the  constellation,  which  puts  high  demands  on  the 
attitude  and  position  control  system.  Darwin  will  be  placed  further  away,  at  a  distance  of  1.5  million 
kilometers  from  Earth,  in  the  opposite  direction  from  the  Sun  (Earth-Sun  Lagrangian  Point  L2  -Fig.8). 


Fig.  7  Darwin  spacecraft  (Artist's  view.  ESA  courtesy) 


Earth 


Lagrange  x  \#; 

points  _  L5 

Sun-Earth  Distance  =  1  AU  =  150,000,000  km 


Fig.  8  Earth-Sun  Lagrangian  Points  and  Darwin  spacecraft  location 

Each  telescope  flyer  is  cylindrically  shaped  (2  m  diameter,  2  m  height)  and  weighs  500  kg.  In  order  to 
protect  the  instrument  from  the  sunlight,  it  is  equipped  with  a  sunshield  modeled  with  6  large  flexible 
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beams  (4  m  long  and  7  kg)  attached  to  the  rigid  structure  (Fig.  9;  beam  end-point  coordinates  in  brackets). 
The  main  mechanical  characteristics  of  the  Darwin-type  Flyer  are  summarized  in  Table  II. 

For  every  beam  (Fig.  9),  two  different  frequencies  for  the  first  modes  along  Y  and  Z  beam  axes  are 
considered.  Their  frequency  can  vary  from  0.05  Hz  to  0.5  Hz,  with  a  nominal  value  of  0.1  Hz,  and  their 
damping  can  vary  from  0.1%  to  1%,  with  a  nominal  value  of  0.5%.  As  regards  spacecraft  mass  and  inertia, 
the  corresponding  uncertainty  around  their  nominal  value  is  of  5%. 

Based  on  the  previous  description  and  using  a  mechanical  modeling  formulation  for  multiple  flexible 
appendages  of  a  rigid  body  spacecraft,  the  open-loop  transfer  function  matrix  representation  of  the 
Darwin-type  Flyer  is  given  in  (50)  and  Fig.  10: 


x(s) 

>11(5)  Pu(s)  Pn(s)  Pu(s)  Pis(s)  P\AS) 

~Fx(s) 

As) 

Pl\(s)  P2l(s)  Pliis)  p24(s)  p25(s)  P26(s) 

Fy(s) 

z(s) 

II 

co 

& 

Co 

X 

II 

P31  (s)  p32  (5)  p33(s)  p34(s)  p3s(s)  p3fr(s) 

Fz(s) 

<p0) 

P4l  (5)  £>42(5)  ^43(5')  P44  (s)  P45(s)  P46(s) 

T<?(s) 

0(5) 

P 51  C5)  P52  (^)  P53  (s)  P54  (5‘)  P 55  (s)  P 56  (5) 

W 

y(s)_ 

_P6l(s)  P62{S)  p63(s )  pM(s)  p65{s)  p66(s)_ 

Tv(s) 

where  x,  y,  z  are  the  position  coordinates;  cp,  0,  \| /  are  the  corresponding  attitude  angles;  Fx,  Fy,  Fz  are  the 
force  inputs;  Ttp,  T0,  T ^  are  the  torque  inputs;  and  where  every  pi](s),  i,  j  =  1,...,6,  is  a  50th  order  Laplace 
transfer  function  with  uncertainty. 


Fig.  9  Darwin  type  6  DOF  satellite  model 


The  Bode  diagram  of  the  plant  (Fig.  10)  shows  the  dynamics  of  the  36  matrix  elements.  Each  of  them  and 
the  MIMO  system  (matrix)  are  minimum  phase.  The  flexible  modes  introduced  by  the  appendages 
(second-order  dipoles)  affect  all  the  elements  around  the  frequencies  oo  =  [0.19,  10]  rad/sec.  The  diagonal 
elements  /?ii(s),  i  =  1,...,6,  and  the  elements  /?i5(s),  /?5i(s),  Pia{s)  and p^iis)  are  mainly  double  integrators 
plus  the  flexible  modes. 
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Table  II.  Mechanical  characteristics  of  the  Darwin-type  Flyer  model 


Parameter 

Value 

Satellite  body  mass 

500  kg 

Cylinder  dimensions 

2  m  diameter,  2  m  height 

Inertia  tensor  of  satellite  in  control  frame  at  satellite 
Centre  of  Mass 
(without  reflector) 

Darwin  body  ]coM  Cont 

'250  0  0  ~ 

0  250  0 

0  0  250_ 

kg-m2 

Inertia  tensor  of  satellite  in  control  frame  at  satellite 
Centre  of  Mass 
(with  reflector) 

[^Darwinbody  ]c0]Vr  1  Cont  — 

"509  0  0  ' 

0  509  0 

0  0  684 

kg-m2 

Position  of  Centre  of  Mass  in  control  frame  at  satellite 
Centre  of  Mass 

[0,  0,  0]  m 

Sunshield  mass 

7  kg  *  6  beams  =  42  kg 

Beam  length 

4  m 

Frequency  Bode-Plant 


From:  ln(1)  From:  ln(2)  From:  ln(3)  From:  ln(4)  From:  ln(5)  From:  ln(6) 


Frequency  (rad /sec) 


Fig.  10  Darwin-type  flyer  dynamics 


The  block  diagram  of  the  control  system  is  shown  in  Fig.  11.  The  sensor  module  represents  both  the  OPD 
(Optical  Pathlength  Differences)  Fringe  Tracker  sensor  and  the  FPM  (Fine  Pointing  Metrology)  sensor, 
which  measure  the  satellite  position  and  attitude,  respectively.  The  actuators,  FEEP  (Field  Emission 
Electric  Propulsion)  thrusters,  are  a  type  of  electrostatic  propulsion  that  provides  very  small  and  precise 
actuation  (Table  III). 
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Fig.  11  General  6x6  satellite  control  loop 

The  external  disturbances  acting  on  the  satellite  (gravity  gradient  and  solar  pressure),  although  very  small, 
are  also  modeled  as  forces  and  torques  along  the  3  axes.  The  gravity  gradient  is  modeled  as  a  constant  bias 
and  the  solar  pressure  is  represented  as  a  white  noise  perturbation  (Table  III). 


Table  III.  Characteristics  of  sensors,  actuators  and  external  disturbances 


Name 

Characteristics 

Values 

Fine  Pointing 

Metrology  (FPM) 

For  very  precise  relative  attitude 
measurements 

0meas  —  ©true  WNfPM 

WNfpm  =  Attitude  white  noise. 

PSD  of  10.66  mas/  a/ Hz  along  the  3  axes 

Attitude  range: [-40;  40]  arcsec 

Sampling  frequency:  1Hz 

Optical  Pathlength 
Differences  (OPD) 
Fringe  Tracker 

For  precise  3 -axis  measure  of  position 
Xmeas  Xtrue  +  WNoPD 

WNopd  -  Position  white  noise. 

PSD  of  2  nm/  a/Hz  along  the  3  axes 

Attitude  range:  [-1;  1]  pm 

Sampling  frequency:  1Hz 

FEEPS  actuators 

For  very  small  and  precise  actuation 

FfEEP  =  F  commanded  +  WNpEEP 

TfeEP  =  Tcommanded  +  WNtFEEP 

Force  model:  WNfeep  =  Force  white  noise. 

PSD  of  0.1  pN  /  a/Hz  along  the  3  axes,  which  can  eventually 

vary  up  to  0.5  pN  /  a/Hz  . 

Force  range:  [-150;  150]  pN 

Torque  model:  WNT  feep  =  Torque  white  noise.  PSD  of  0.1 
pNm  /  a/ Hz  along  the  3  axes,  which  can  eventually  vary  up 

to  0.5  pNm  /  a/Hz  . 

Torque  range:  [-150;  150]  pNm 

Sampling  frequency:  1Hz 

Solar  Pressure 

Fsun  Bp_Sun  W^NF  Sun 

Tsun  —  B'r  Suti  NTSun 

Force  model:  WNf  sun  =  Force  white  noise. 

PSD  of  0.05  pN  /  a/Hz  along  the  3  axes 

Bf  Sun  =10  pN 

Torque  model:  WNT  sun  =  Torque  white  noise.  PSD  of  0.1  pNm 

/  a/ Hz  along  the  3  axes 

Bt  Sun  =  10  pNm 

Gravity  Gradient 

FGrav 

TGrav 

Force  model:  FGrav  =  0.03  pN  along  the  3  axes 

Torque  model:  TGrav  =  0.03  pNm  along  the  3  axes 

0  =  attitude.  X  =  position.  F  =  Force.  T  =  Torque 


The  original  dynamics  benchmark  simulator,  provided  by  ESA  and  implemented  under  Matlab/Simulink, 
integrates  all  those  elements  constituting  the  whole  satellite  control  system:  sensors,  actuators,  dynamics, 
disturbances,  etc.  (Fig.  11).  For  each  performance  evaluation  campaign,  300  random  dynamics  within  the 
uncertainty  (Monte-Carlo  analysis)  are  generated  to  evaluate  the  performance  of  the  controllers. 
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Table  IV.  Darwin-type  Flyer  requirements 


Objective 

Numerical  Requirement 

Astronomical  Requirements 

Position  accuracy 

Maximum  absolute  value: 

1  |im  for  all  axes 

Standard  deviation: 

0.33  pm  for  all  axes 

Pointing  accuracy 

Maximum  absolute  value: 

25.5  mas  for  all  axes  (3  a) 

Standard  deviation: 

8.5  mas  for  all  axes  (la) 

Engineering  Requirements 

Bandwidth 

~  0.01  Hz  for  all  axes 

Saturation  limits 

Maximum  force:  150  pN 
Maximum  torque:  150  pNm 

Rejection  of  high  frequency  noises 
(from  measurement  and  actuation) 

High  roll-off  after  the  bandwidth 

Control  Requirements 

Stability  margins 

max\T(jco) 1  <  2 

CO 

max\s(jco^  <  2 

Loop  interaction 

Minimum 

Rejection  of  flexible  modes 

Maximum 

Controller  complexity  and  order 

Minimum 

3.2  Control  objectives 

The  main  objective  of  the  spacecraft  is  to  fulfill  some  astronomical  requirements  that  demand  to  keep  the 
flying  telescope  pointing  at  both  the  observed  space  target  and  the  central  hub-satellite.  This  set  of 
specifications  leads  to  some  additional  engineering  requirements  (bandwidth,  saturation  limits,  noise 
rejection,  etc.)  and  also  needs  some  complementary  control  requirements  (stability,  low  loop  interaction, 
low  controller  complexity  and  order,  etc.)  -Table  IV-. 

3.3  Non-diagonal  MIMO  QFT  Controller  Design 

The  sequential  non-diagonal  MIMO  QFT  methodology  previously  described  in  Section  2  [2-7]  is  applied 
here  to  control  the  position  and  attitude  of  the  Darwin-type  Flyer. 


3.3.1  Relative  Gain  Array  Interaction  Analysis  -Step  A- 

The  Relative  Gain  Array  (RGA)  of  a  non-singular  square  matrix  P  is  a  scale-invariant  measure  of 
interactions.  Its  original  definition  introduced  by  Bristol  [11]  was  only  proposed  for  steady  state  (oe>  =  0 
rad/sec).  However,  the  RGA  can  also  be  computed  frequency-by-frequency  (51)  and  used  to  assess  the 
interaction  at  frequencies  other  than  zero  [10].  In  most  cases,  the  value  of  RGA  at  frequencies  close  to 
crossover  is  the  most  important  one. 


RGA(]co)=  [^y (j &>)]  =  -PQ ffl) ® (p  ](]oj))T 


(51) 


where  ®  denotes  element-by-element  multiplication  (Schur  product).  Further  information  on  how  to 
interpret  the  RGA  results  and  select  pairings  can  be  found  at  [10,  11]. 

The  6x6  (position  and  attitude)  dynamic  model  of  the  Darwin-type  spacecraft  with  large  flimsy 
appendages  has  been  analyzed  by  using  the  RGA  method  as  a  function  of  frequency  and  for  the  whole  set 
of  parameter  uncertainty.  Although  the  matrix  obtained  by  means  of  (51)  is  a  complex  one,  only  the 
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absolute  values  are  presented.  By  examining  the  corresponding  matrices  at  each  frequency,  it  is  observed 
that  the  steady  state  results  extend  through  low  frequency  up  to  0.19  rad/sec.  As  a  representative  example 
within  this  range,  (52)  shows  the  results  for  the  most  coupled  plant  within  the  uncertainty  at  00  =  6.28- 10'4 
rad/sec.  According  to  it,  the  pairing  should  be  done  through  the  main  diagonal  of  the  matrix,  which 
contains  positive  RGA  elements,  and  the  elements  g\s(s%  gi4(s),  g4i(s),  gsi(s)  should  also  be  considered 
relevant. 


RGA 


0=6.28-1 0“4  rad/sec 


1.0064  0  0  0  0.0064  0 

0  1.0064  0  0.0064  0  0 

0  0  1  0  0  0 

0  0.0064  0  1.0064  0  0 

0.0064  0  0  0  1.0064  0 

0  0  0  0  0  1 


(52) 


If  the  analysis  is  performed  at  high  frequency,  it  produces  the  same  concluding  results  in  the  entire 
spectrum  starting  at  3.51  rad/sec. 

So  far,  the  retained  compensator  elements  would  be  those  of  the  RGA  matrix  marked  in  bold  in  (52). 
Nevertheless,  as  aforementioned,  the  RGA  elements  increase  and  more  interactions  are  founded  in  the 
interval  of  frequencies  where  the  flexible  modes  of  the  satellite  mostly  affect  (i.e.  00  =  [0.19-3.51]  rad/sec), 
as  can  be  seen  in  (53)  and  (54)  for  the  most  coupled  plants  at  co  =  0.8  rad/sec  and  1  rad/sec,  respectively. 


^^4  (<^=0.8  rad/sec) 


(a>=\  rad/sec) 


1.1674 

0.0001 

0.0001 

1.0180 

0.0000 

0.0000 

0.0012 

0.0305 

0.3720 

0.0013 

0.0001 

0.0004 

0.9899 

0.0002 

0.0002 

0.9082 

0.0000 

0.0000 

0.0001 

0.0050 

0.0003 

0.0009 

0.0102 

0.0963 

0.0000 

0.0012 

0.0000 

0.0305 

4.8592 

0.3468 

0.3468 

3.2030 

4.1456 

2.3865 

0.0000 

0.0006 

0.0000 

0.0001 

0.0000 

0.0050 

3.1307 

2.3373 

2.3373 

12.5674 

0.1110 

9.4613 

0.0000 

0.0690 

0.3720  0.0001 
0.0013  0.0004 
4.1456  0.0000 
2.3865  0.0006 
7.2470  0.0008 
.0008  1.0009 


0.0003 

0.0102 

0.0009 

0.0963 

0.1110 

0.0000 

9.4613 

0.0690 

10.2042 

0.0302 

0.0302 

0.7970 

(53) 


(54) 


3.3.2  Controller  Structure 


In  accordance  with  the  above  RGA  results  and  taking  into  account  the  requirement  of  minimum  controller 
complexity  and  order  (Table  IV),  a  first  compensator  structure  consisting  of  six  diagonal  elements  and 
four  off-diagonal  elements  is  chosen  as  the  most  suitable  one  (55). 


g\\(s)  0  0  0  £15  0)  0 

0  £22  C^)  0  £24  0)  0  0 

0  0  £33^)  0  0  0 

0  §42  0s)  0  §44  (s)  0  0 

§5l(J)  0  0  0  g55(.s)  0 

0  0  0  0  0  g66(s) 


(55) 
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From  this,  four  independent  compensator  design  problems  have  been  adopted,  two  SISO  and  two  2x2 
MIMO  problems:  [g33(s)]  and  [geeO)];  [gn(s)  gisO)  ;  gsi(s)  gssO)]  and  [g22(s)  g24(s)  ;  g4i(,s)  g44^)], 
respectively.  The  SISO  problems  will  be  considered  from  the  classical  SISO  QFT  point  of  view,  while  the 
two  2x2  MIMO  subsystems  will  be  studied  through  the  non-diagonal  MIMO  QFT  methodology.  The 
coupling  detected  in  the  range  of  frequencies  of  the  flexible  modes  will  be  considered  in  the  course  of  the 
design  procedure  through  more  demanding  specifications  with  respect  to  disturbance  rejection.  Provided 
the  performance  results  were  not  satisfactory,  then  the  compensator  structure  should  be  filled  up  with 
additional  off-diagonal  compensators  consistent  with  (53)  and  (54):  g34,  g35,  g43,  g45,  g53  andg54  elements. 

3.3.3  Robust  Closed-Loop  Specifications 

The  applied  non-diagonal  and  diagonal  MIMO  QFT  techniques  design  each  loop  individually,  including 
the  multivariable  characteristic  by  means  of  their  corresponding  equivalent  plant.  So,  the  robust  closed- 
loop  specifications  are  defined  in  terms  of  SISO  expressions  for  both  SISO  and  MIMO  subsystems. 

Since  these  methodologies  are  frequency  domain  techniques,  there  is  obviously  a  need  for  translating  time 
domain  requirements  (Table  IV)  into  the  frequency  domain.  The  original  specifications  in  Table  IV 
impose  maximum  and  standard  deviation  values  on  the  position  and  attitude  time  error  signals,  as  well  as 
actuator  forces  and  torques.  Their  translation  into  the  frequency  domain  is  based  on  control-ratio  models 
[9],  and  takes  into  account  the  expected  external  disturbances  on  the  Darwin-type  flyer,  the  spacecraft 
flexible  modes  and  the  coupling  among  loops.  As  a  result,  four  Type  of  specifications  are  defined  to 
calculate  the  QFT  bounds:  Type  1:  Robust  stability;  Type  2:  Robust  sensitivity;  Type  3:  Robust 
disturbance  rejection  at  plant  input;  and  Type  4:  Robust  control  effort  attenuation. 

The  notation  used  for  the  signals  in  the  following  description  of  specifications  refers  to  the  scheme  of  the 
generic  MIMO  subsystem  presented  in  Fig.  12.  The  compensators  have  been  designed  within  the  set  of 
frequencies  of  interest  co  =  [6.28- 10"4,  62.8]  rad/sec. 


Fig.  12  Structure  of  a  2  Degree  of  Freedom  MIMO  System 
Type  1:  Robust  Stability  specification 

This  specification,  shown  in  (56),  is  stated  to  guarantee  a  robust  stable  control.  All  the  required  values, 
displayed  loop  by  loop  in  (57)  and  (58),  imply  at  least  1.54  (3.75  dB)  gain  margin  and  at  least  49.25° 
phase  margin.  The  specification  corresponds  not  only  to  the  closed-loop  transfer  function  y^/r^s),  but 
also  to  transfer  functions  yi(s)/?7i(s)  and  Wi(s)/vi(s).  Hence  this  condition  additionally  imposes  the 
requirements  on  sensor  noise  attenuation,  disturbance  rejection  at  plant  input  and  flexible  modes. 


*q(  d 

-i  ( \ 

[Pn  wj, 

i  8h{s) 

1+[/VeCs\ 

Jr1^) 

<  Sx{(o) 


(56) 


where  [pii*e(^)]i_1  is  the  inverse  of  the  equivalent  plant  (45),  which  corresponds  to  /?ii(s)  in  SISO  designs. 
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Loops  1,  2  and  3: 


Loops  4,  5  and  6: 


Si(co)  =  1.85  ;  Vco 


Si(a>)= 


0.1687 


sz  +  0.4s +  0.0912 


Vco 


(57) 

(58) 


Type  2:  Sensitivity  reduction 

The  main  objective  of  this  specification,  (59)  and  (60),  is  sensor  noise  attenuation  and  reduction  of  the 
effect  of  the  parameter  uncertainty  on  the  closed-loop  transfer  function.  It  corresponds  to  ei(s)/n[(s)  and 
[dAi(s,)/Ai(s')]  /  [6pn(s)/pn(s)]  transfer  functions. 


1 


1  + 


<  S2(co) 


(59) 


All  loops:  S2(co)  =  2  ;  Vco 


(60) 


Type  3:  Rejection  of  disturbances  at  piant  input 

Solar  pressure  perturbation  and  gravity  gradient  are  considered  to  affect  at  plant  input  in  the  form  of  both 
force  and  torque.  The  purpose  of  this  specification  (61),  which  corresponds  to  e\(s)/v\(s)  and  y\(s)/v;(s) 
transfer  functions,  is  to  attenuate  the  effect  of  plant  input  disturbances  on  the  control  error  and  the  output 
signal.  Thus,  a  high  gain  is  required  in  the  low  frequency  band,  (62)  to  (64).  Besides,  since  Vi(s)  also 
represents  the  flexible  modes,  special  attention  is  paid  to  their  frequency  range  mainly  to  accomplish  the 
attitude  requirements. 


+ 

’ii*e(5)]i”+ii(5) 

<  83{cd) 


(61) 


Loop  3: 
Loops  4,  5  and  6: 


Loops  land  2:  S3(co)= 
S3(co)= 


0.21553  (s  + 0.385) 

(s  +  0.307)  (s  +  6. 18)  (s2  +  0.4i  +  0.0912) 


;  Vco 


0.313  (s-0. 01705)  (s2  +0.009974s  +  5.104-10~5) 


(s-0.01813)  (s2  +  0.02554s +  0.0004754) 


(s  +  0.2)  (s  +  0.186)  (s  +  0.2044)  (s  +  0.003892)  (s2  +  0.06014s  +  0.02736) 
(s  +  0.007333)  (s  +  0.445)  (s2  +  0.07904s  +  0.00326)(,s2  +  0.2352s  +  0.098l) 


Vco 


(62) 


(63) 


;  Vco  (64) 


Type  4:  Control  signal 


Because  of  saturation  limits,  control  signal  movements  should  be  kept  reasonably  small  despite 
disturbances  coming  from  actuators  and  sensors.  This  specification,  (65),  corresponds  to  Ui{s)lnis)  transfer 
function  and  is  depicted  in  (66)-(68). 


1  + 


Sn(S) 


<  S4(co) 


(65) 
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Loops  1  and  2:  S4(co)= 


557.1  (s  +  5) 


51 


Vco 


Loop  3: 


S4(a>)  = 


s'  +  3.23s +  6.5 
106.9210  (s  +  0.55)(^2  +  0.04s +  0.13) 


Loops  4,  5  and  6:  S4(a>)= 


(s  +  1.4)2  (^2  +  0.1227s +  0.097) 

4.026  (s2  - 0. 1854s  +  0.203)(.s2  +  0.04s  +  0.504) 


;  Vco 


(s2  +  0.305s  +  0.056)(s2  +0.1 15s  +  0.095) 


;  Vco 


(66) 


(67) 


(68) 


Reducing  coupling  effects  as  much  as  possible 

The  coupling  effects  from  other  axes  can  be  considered  as  part  of  the  disturbances  acting  at  the  input  of 
the  equivalent  SISO  plant.  The  way  of  designing  the  non-diagonal  elements  of  the  matrix  compensator 
deals  with  the  aim  of  minimizing  the  off-diagonal  elements  of  the  coupling  matrix  (32). 


3.3.4  SISO  Design  Problems:  £33(3),  g66(s) 

Compensators  g33(s)  and  g66(s)  are  independently  designed  by  using  classical  SISO  QFT  [1]  to  satisfy  the 
robust  stability  and  robust  performance  specifications  stated  in  Section  3.3.3  for  every  plant  within  the  set 
of  uncertain  plants.  The  corresponding  QFT  bounds  and  the  nominal  case  of  the  designed  open-loop 
transfer  functions  Zii(s)  =  /?ii(s)  gn(s),  i  =  3,  6,  are  plotted  on  the  Nichols  Chart  for  some  of  the  most 
relevant  frequencies  in  Fig.  13(a)  and  13(b)  for  loops  3  and  6  respectively.  Both  designs  satisfy  not  only 
their  respective  bounds  but  also  the  Nyquist  encirclement  condition,  and  no  RHP  pole-zero  cancellations 
occur  between  g33(s)  and  ^>33(5),  nor  between  g66(s)  and  pee(s).  The  Bode  plot  of  each  compensator  can  be 
found  in  Section  3.6,  where  g33(s)  [Fig.  18(a)]  and  g66(s)  [Fig.  18(b)]  are  represented  in  solid  line  in 
comparison  with  the  H-infinity  design  (dashed  line)  introduced  in  Section  3.5.  The  QFT  compensator 
expressions  are  presented  in  Section  3.5. 


Loop-Shaping  Loop- Shaping 


(a)  (b) 


Fig.  13  Loop-shaping  (a)  L33(s)  =  p33(s)  g33(s).  (b)  L66(s)  =  p66(s)  g66(s) 
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3.3.4.1  First  MIMO  Problem:  gn(s),  £51 0),  gss(s),  gis(s)  Design 

The  compensator  for  this  2x2  MIMO  problem  has  been  designed  by  applying  the  non-diagonal  MIMO 
QFT  methodology  developed  by  Garcia-Sanz  et  al.  [1-7]  and  outlined  in  Section  2.  In  this  particular  case, 
the  plant  to  be  controlled  is  composed  of  the  following  elements  coming  from  the  original  6x6  Darwin- 
type  spacecraft  model  Pl5(s)  =  [pn(V)  ;  /?5i(X)  P5s(s)],  whose  inverse  matrix  is  P15*^)  =  [P15^)]  1  = 

bn*0)  pi5\s)  ;p5i\s)  p55*(s)]. 

Step  A:  Arrangement  of  the  system 

First,  the  methodology  adopts  the  structure  and  the  pairing  of  inputs  and  outputs  given  by  the  RGA 
technique  in  (55)  and  arranges  the  plant  inverse  matrix  P15*(s )  so  that  the  inverse  of  the  first  diagonal 
element  in  this  matrix  has  the  smallest  phase  margin  frequency  [1].  In  some  cases,  arbitrarily  picking  the 
wrong  order  of  the  loops  could  lead  to  the  non-existence  of  a  solution.  In  the  present  problem,  the 
bandwidth  of  the  loops  is  quite  similar.  Then,  any  order  can  be  selected  to  design  the  non-diagonal  MIMO 
QFT  compensators. 

Step  B1:  Design  of  the  diagonal  compensator  g^fs) 

The  diagonal  compensator  gn(s)  is  designed  through  standard  QFT  loop-shaping  [1]  for  the  inverse  of  the 
equivalent  plant  \pn  e(s)]i  =  pn  (s')  to  fulfill  the  robust  stability  and  robust  performance  specifications 
determined  in  Section  3.3.3  for  every  plant  within  the  set  of  uncertain  plants.  Fig.  14(a)  shows  the  nominal 
case  of  the  designed  open-loop  transfer  function  Zn^)  =  [pii^CQlf 1  gn(s)  in  bold  solid  line,  which 
satisfies  the  QFT  bounds,  also  represented  in  the  figure.  Additionally,  the  design  fulfils  the  first  two 
sufficient  stability  conditions  (c.l)  and  (c.2)  (Section  2.6).  That  is,  Zn(5)  =  \pn*G(s)]fl  gn(s)  satisfies  the 
Nyquist  encirclement  condition  and  no  RHP  pole-zero  cancellations  occur  between  gn(s)  and  [pii*e(^)]i_1. 
The  Bode  plot  for  the  obtained  compensator  gn(s)  is  presented  in  Fig.  19(a)  (solid  line)  together  with  the 
design  of  the  H-infinity  approach. 


Loop-Shaping  Loop-Shaping 


Fig.  14  Loop-shaping  (a)  Ln(s)  =  [pn^sjlr1  gn(s).  (b)  L55(s)  =  [p55*e(s)]2_1  gss(s) 

Step  Cp  Design  of  the  non-diagonal  compensator  g51(s) 

The  non-diagonal  compensator  gs\(s)  is  designed  to  minimize  the  (5,1)  element  of  the  coupling  matrix  in 
the  case  of  disturbance  rejection  at  plant  input  (32),  which  gives  the  following  expression: 


RTO-EN-SCI-195 


1  -23 


Beyond  the  Classical  Performance  Limitations  Controlling  Uncertain  MIMO  Systems 


ORGANIZATION 


g5Pl(s)  =  -P5lN(s)  (69> 

where  N  denotes  the  plant  that  minimizes  the  maximum  of  the  non-parametric  uncertainty  radii 
comprising  the  plant  templates  on  the  Nichols  Chart.  Due  to  the  uncertainty,  the  expression  [~p5 1  (s)] 
determines  a  region  in  the  magnitude  and  phase  plots,  where  the  compensator  g5i(s)  is  shaped  following 
the  mean  value  at  every  frequency  oo  e  [0,  0.1]  rad/sec  [see  Fig.  15  with  g5i(s)  interpolating  the  plot].  The 
compensator  Bode  plot  is  compared  in  Fig.  19(c)  with  that  of  the  (5,1)  element  of  the  H-inflnity 
compensator  introduced  in  Section  3.5. 


Magnitude  plot  of  g^($) 


10-3  10 2  10 1  10°  101  102 
Frequency  (red/eec) 


Fig.  15  Magnitude  plot  of  [-psi  (s)]  with  uncertainty  and  gsi(s)  -bold  solid  line- 
Step  B2:  Design  of  the  diagonal  compensator  g5s(s) 

As  in  step  B u  the  diagonal  compensator  g55(s)  is  designed  to  control  the  inverse  of  the  equivalent  plant, 
\P55  e(s)]2  \  which  takes  the  compensator  previously  designed  into  account  (45). 


[/J55^)]2=[/4e(^)]l 


( te 

)1  +b51(s)]1)( 

:te(4) 

l 

p'M] 

1,+kM], 

(70) 


On  the  basis  of  the  robust  specifications  defined  in  Section  3.3.3  for  [p55*e(^)]2  \  and  also  taking  into 
account  the  plant  uncertainty,  the  QFT  bounds  are  computed.  Then,  the  compensator  is  designed  by 
classical  loop-shaping  on  the  Nichols  Chart,  as  is  shown  in  Fig.  14(b).  Not  only  does  the  design  fulfil  the 
bounds  but  also  the  first  two  stability  conditions  of  (c.l)  and  (c.2)  from  Section  2.6.  In  other  words,  L55(s ) 
=  [p55  e(s)]2_1  g55(s)  satisfies  the  Nyquist  encirclement  condition  and  no  RHP  pole-zero  cancellations  occur 
between  g55(s)  and  [p55*e(s)]2  The  Bode  plot  of  g55(s)  is  presented  in  Fig.  19(d). 

Step  C2:  Design  of  the  non-diagonal  compensator  g15(s) 

Due  to  the  requirement  of  minimum  controller  complexity  and  order  (Table  IV),  the  non-diagonal 
compensator  gi5(s)  has  been  set  to  zero.  Anyway,  the  equivalent  expression  to  the  one  used  in  (69), 
gi5°pt(s)  =  -Pi5*N(s)?  could  be  applied  to  cancel  interaction  in  both  directions  in  the  MIMO  subsystem. 

At  this  point,  once  the  whole  controller  of  the  MIMO  subsystem  has  been  determined,  the  last  two 
stability  conditions  mentioned  in  Section  2.6,  (c.3)  and  (c.4),  are  checked.  The  system  is  stable.  Finally, 
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the  non-existence  of  RHP  transmission  zeros  of  P(s )  G(s )  is  checked  by  using  the  Smith-McMillan  form 
over  the  set  of  possible  plants  3P  due  to  uncertainty  [5].  The  non-diagonal  MIMO  QFT  compensator 
expressions  are  presented  in  Section  3.6. 

3.3.4.2  Second  MIMO  Problem:  g 22(2),  g42 (s),  g44 (s),  g24(s)  Design 

The  second  MIMO  problem  consists  of  the  following  elements:  P24^)  =  [puis)  P2a(s)  ;  Pai{s)  /?44(s)]. 
From  the  2x2  plant  inverse  matrix  P24*^)  =  [P24^)]"1  =  \p2i(s)  P2a(s)  ;  p4 2*(s)  p44*(s)]  and  taking  into 
account  the  robust  stability  and  robust  performance  specifications  (Section  3.3.3),  the  non-diagonal  MIMO 
QFT  methodology  is  equivalently  performed  by  following  the  steps  detailed  in  Section  2. 

The  loop-shaping  for  the  diagonal  compensator  elements  g22(s)  and  g44(s)  are  shown  in  Fig.  16(a)  and 
16(b),  respectively.  The  Bode  plots  for  the  four  compensators  are  shown  in  Fig.  20(a),  (b),  (c)  and  (d)  for 
g22(s),  g24(s),  g4i(s),  g44(s),  respectively.  The  2x2  MIMO  subsystem  is  found  to  be  stable  according  to  the 
sufficient  stability  conditions  (Section  2.6).  Finally,  it  is  also  checked  that  no  additional  RHP  zeros  have 
been  introduced  by  the  compensator  [5].  The  non-diagonal  MIMO  QFT  compensator  expressions  are 
presented  in  Section  3.6. 


Loop-Shaping  Loop-Shaping  L, 


(a)  (b) 

Fig.  16  Loop-shaping  (a)  L22(s)  =  [p22*e(s)]r1  g22(s).  (b)  L44(s)  =  [p44*e(s)]2_1  g44(s) 


3.4  Diagonal  MIMO  QFT  Controller  Design 

For  the  sake  of  comparison,  the  sequential  diagonal  MIMO  QFT  methodology  developed  by  Horowitz 
[12]  is  also  applied  to  control  the  position  and  attitude  of  the  Darwin-type  Flyer.  Based  on  the  same  robust 
closed-loop  specifications  defined  in  Section  3.3.3,  this  technique  uses  a  sequential  procedure  similar  to 
the  one  detailed  in  Section  2.5  (Step  B),  where  the  recursive  expression  of  the  equivalent  plant  is  a 
simplified  case  of  (45),  with  gy  (5)  =  0  (i  ^  j). 

For  the  Darwin-type  Flyer,  the  loop-shaping  step  of  the  diagonal  method  requires  the  same  diagonal 
compensators  ga  (5)  as  the  non-diagonal  one.  This  happens  because,  in  this  case,  in  the  middle  and  high 
frequency  range  the  off-diagonal  elements  gij(s)  (i  ±  j)  of  the  non-diagonal  controller  have  less  relative 
influence  than  the  corresponding  pX]  (s)  elements  in  the  equivalent  plant  (45).  Differences  between  both 
MIMO  QFT  controllers  arise  in  the  low  frequency  range,  as  can  be  observed  in  Fig.  17.  The  C2(5,l) 
element  of  the  coupling  matrix  for  disturbances  at  plant  input  (32)  is  plotted  for  a  representative  plant  case 
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and  for  the  three  controllers:  non-diagonal  and  diagonal  MIMO  QFT  and  H-infinity  designs.  For  the 
frequency  range  cd  e  [0,  0.1]  rad/sec  it  is  shown  that  C2(5,l)non-diag  qft  <  C2(5,l)diag  qft  <  C2(5,l)H-mfmity, 
which  explains  why  the  non-diagonal  MIMO  QFT  improves  the  diagonal  MIMO  QFT  and  the  H-infinity 
controller  results  under  low  frequency  external  disturbances  (see  Section  3.6.2). 


Coupling  matrix  for  disturbance  rejeclion  at  plant  input:  C2l*j/I) 


Fig.  17  Element  (5,1)  of  the  coupling  matrix  C2'.  non-diagonal  MIMO  QFT  in  solid  line,  diagonal  MIMO  QFT  in 

dotted  line,  H-infinity  in  dashed  line 


3.5  Controllers 

The  notation  adopted  for  transfer  function  expressions  denotes  the  steady  state  gain  as  a  constant  without 
parenthesis;  simple  poles  and  zeros  as  (co),  which  corresponds  to  (s/co  +1)  denominator  and  numerator, 
respectively;  poles  and  zeros  at  the  origin  as  (0);  conjugate  poles  and  zeros  as  [£,  ;  oon],  with  ((s/con)2  + 
(2^/cOn)  5+1)  denominator  and  numerator,  each;  ^-multiplicity  of  poles  and  zeros  as  an  exponent  ( )”. 


The  non-diagonal  MIMO  QFT  compensator  consists  of  the  following  eight  elements:  gn(s)  =  g22(s)  = 
{31.5  (0.6194)  (0.2138)  (0.1663)  (0.1649)}  /  {(0.666)  (0.4982)  (0.07526)  [0.676;  1.479]}  ;  g5l(s)  =  - 
g4i(s )  =  {42.4  (0)2  }  /  {(0.3)3};  g15(s)  =  g24(s)  =  0;  g33(s)  =  {125  (0.13)  (0.057)  [0.07019;  0.3565]  [1; 
0.02]}  /  {(1.48)  (0.7875)  (0.2)  (0.004)  (0.00246)  [0.18;  0.314]};  g44(s)  =  {2.242  (0.03412)  [0.08644; 
0.7114]  [0.1131;  0.3414]  [0.1145;  0.2604]  [0.008792;  0.2593]  [0.7;  0.0052]  [1;  0.0007]}  /  {(0.9776)  (0.8) 
(0.0005)2  [0.2451;  0.2708]  [0.297;  0.2673]  [-0.0005;  0.254]  [0.14;  0.252]  [0.7;  0.0045]};  g55(s)  =  {2.242 
(0.13)  (0.03)  [0.1079;  0.7099]  [0.07069;  0.341]  [0.03;  0.2593]  [0.7;  0.0052]  [1;  0.0007]}  /  {(0.9776)  (0.8) 
(0.12)  (0.0005)2  [0.2451;  0.2708]  [-0.0008;  0.254]  [0.3;  0.25]  [0.7;  0.0045]};  g66(s)  =  {2.242  (0.02584) 
[0.08644;  0.7114]  [0.1131;  0.3414]  [0.1145;  0.2604]  [0.008792;  0.2593]  [0.7;  0.0052]  [1;  0.0007]}  / 
{(0.9776)  (0.8)  (0.0005)2  [0.2451;  0.2708]  [0.297;  0.2673]  [-0.0007;  0.254]  [0.1687;  0.241]  [0.7; 
0.0045]}. 


The  diagonal  MIMO  QFT  compensator  consists  of  the  same  diagonal  elements  ga(s)  as  the  non-diagonal 
compensator  abovementioned,  and  gij(s)  =  0,  i  +  j. 

The  main  elements  of  the  1-DOF  H-infinity  compensator  are  shown  in  Figs.  18-20.  Their  dc  gains  stay 
within  the  range  [-15  dB,  26  dB].  The  remaining  26  elements  present  a  very  low  gain,  going  from  -260  dB 
to  -330  dB. 
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Darwin  Controllers  Comparison  :  G(33) 


Darwin  Controllers  Comparison  :  Gi]Sj6) 


(a) 


(b) 


Fig.  18  Bode  Diagram  Compensators:  non-diagonal  and  diagonal  MIMO  QFT  in  solid  line, 
H-infinity  in  dashed  line,  (a)  g33(s),  (b)  c/eefs) 


Darwin  Controllers  Comparison  :  G(1,1) 


Frequency  (rad/sec) 

(a) 


Darwin  Controllers  Comparison  :  G(1 .5) 


(b) 


(c) 


(d) 


Fig.  19  Bode  Diagram  Compensators:  non-diagonal  MIMO  QFT  in  solid  line  [also  diagonal  MIMO  QFT  for  gn(s)  and  g55(s)], 

H-infinity  in  dashed  line,  (a)  gn(s),  (b)  gi5(s),  (c)  g5i(s),  (d)  g55(s) 
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Frequency  (rsd/sec)  Frequency  (raCWsec) 

(C)  (d) 


Fig.  20  Bode  Diagram  Compensators:  non-diagonal  MIMO  QFT  in  solid  line  [also  diagonal  MIMO  QFT  for 
g22<s)  and  g44(s)],  H-infinity  in  dashed  line,  (a)  g22(s),  (b)  g24(s),  (c)  g42(s),  (d)  g44<s) 


3.6  Comparative  evaluation 

This  section  shows  a  comparative  analysis  of  the  sequential  non-diagonal  MIMO  QFT  controller,  designed 
above  for  the  6x6  Darwin-type  Flyer,  with  both  sequential  diagonal  MIMO  QFT  and  H-infinity 
controllers.  First,  comparative  Bode  plots  of  the  compensators  are  shown.  Then,  time  performance  results 
are  presented  (astronomical  requirements),  followed  by  open-loop  bandwidth,  and  forces  and  torques 
comparison  (engineering  requirements).  Finally,  the  stability  objectives  and  the  order  of  each  compensator 
are  analyzed  (control  requirements). 

3.6.1  Compensators  Bode  Plots 

The  Bode  plots  are  presented  for  the  compensators  of  the  non-diagonal  MIMO  QFT  (solid  line)  in 
comparison  with  those  of  the  H-infinity  (dashed  line).  Note  that,  in  this  case,  the  diagonal  MIMO  QFT 
method  yields  the  same  diagonal  compensators  as  the  non-diagonal  MIMO  QFT  technique.  Fig.  18 
presents  the  results  for  the  two  SISO  subsystems  g33(s)  and  g6 6(s),  (a)  and  (b)  respectively.  Fig.  19  plots 
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the  compensators  of  the  2x2  MIMO  subsystem  composed  of  gn(s),  gi5(s),  gsi(s)  and  g55(s)  elements.  The 
g22(^),  g24(^),  g42(s)  and  g44(s)  compensator  elements  that  conform  the  other  2x2  MIMO  subsystem  are 
shown  in  Fig.  20.  Note  that  gi5(s)  and  g24(s)  have  been  set  to  zero  in  the  non-diagonal  MIMO  QFT  design. 
Additionally,  according  to  the  RGA  results  in  (55)  the  remaining  elements  of  the  controller  matrix  G(s ) 
designed  with  this  technique  equal  zero.  By  contrast,  those  26  elements  present  a  non-zero,  although  very 
small,  magnitude  response  when  they  are  designed  with  the  H-inflnity  technique.  Finally,  the  off-diagonal 
elements  of  the  diagonal  MIMO  QFT  are  obviously  zero. 


Table  V.  Time  simulation  performance  with  the  three  controllers 


Specification 

Requirement 

Benchmark 

Non-diagonal 
MIMO  QFT 

Diagonal 
MIMO  QFT 

H-infinity 

Controller 

Controller 

Controller 

1 

Maximum  Position 

<  1  pm 

B1 

0.0131 

0.0131 

0.0293 

Error  X  (pm) 

B2 

0.0816 

0.0816 

0.511 

Maximum  Position 

<  1  pm 

B1 

0.0120 

0.0120 

0.0299 

z 

Error  Y  (pm) 

B2 

0.0120 

0.0120 

0.0299 

T 

Maximum  Position 

<  1  pm 

B1 

0.0288 

0.0288 

0.0292 

J 

Error  Z  (pm) 

B2 

0.0288 

0.0288 

0.0292 

A 

Maximum  Attitude 

<  25.5  mas 

B1 

25.27 

25.31 

25.95 

4 

Error  X  (mas) 

B2 

25.27 

25.31 

25.95 

c 

Maximum  Attitude 

<  25.5  mas 

B1 

22.91 

22.99 

23.21 

J 

Error  Y  (mas) 

B2 

22.55 

23.75 

28.91 

Maximum  Attitude 

<  25.5  mas 

B1 

21.15 

21.15 

22.84 

O 

Error  Z  (mas) 

B2 

21.15 

21.15 

22.84 

7 

Std.  Deviation  of 
Position  Error  X 
(lim) 

<  0.33  pm 

B1 

B2 

0.00275 

0.0511 

0.00276 

0.0511 

0.00686 

0.341 

8 

Std.  Deviation  of 
Position  Error  Y 
(nm) 

<  0.33  pm 

B1 

B2 

0.00265 

0.00265 

0.00266 

0.00266 

0.00722 

0.00722 

9 

Std.  Deviation  of 

<  0.33  pm 

B1 

0.00668 

0.00668 

0.00691 

Position  Error  Z  (pm) 

B2 

0.00668 

0.00668 

0.00691 

10 

Std.  Deviation  of 
Attitude  Error  X 
(mas) 

<  8.5  mas 

B1 

B2 

5.57 

5.57 

5.57 

5.57 

5.68 

5.68 

Std.  Deviation  of 

B1 

B2 

5.76 

5.80 

5.76 

5.85 

6.01 

8.23 

11 

Attitude  Error  Y 
(mas) 

<  8.5  mas 

Std.  Deviation  of 

B1 

B2 

4.83 

4.83 

4.83 

4.83 

5.00 

5.00 

12 

Attitude  Error  Z 
(mas) 

<  8.5  mas 

13 

Maximum  Actuator 
Force  Command  X 
(N) 

<  1.5e-4  N 

B1 

B2 

1.94e-6 

3.94e-6 

1.94e-6 

3.94e-6 

7.42e-7 

3.31e-6 

14 

Maximum  Actuator 
Force  Command  Y 
(N) 

<  1.5e-4  N 

B1 

B2 

1.86e-6 

1.86e-6 

1.86e-6 

1.86e-6 

6.68e-7 

6.68e-7 

15 

Maximum  Actuator 
Force  Command  Z 
(N) 

<  1.5e-4  N 

B1 

B2 

5.94e-7 

5.94e-7 

5.94e-7 

5.94e-7 

5.61e-7 

5.61e-7 

Maximum  Actuator 

B1 

B2 

8.68e-7 

8.68e-7 

8.71e-7 

8.71e-7 

1.03e-6 

1.03e-6 

16 

Torque  Command  X 
(Nm) 

<  1.5e-4  N  m 

Maximum  Actuator 

B1 

B2 

1.05e-6 

1.06e-6 

1.05e-6 

1.06e-6 

1.15e-6 

1.16e-6 

17 

Torque  Command  Y 
(Nm) 

<  1.5e-4  N  m 

18 

Maximum  Actuator 
Torque  Command  Z 
(Nm) 

<  1.5e-4  N  m 

B1 

B2 

1.08e-6 

1.08e-6 

1.08e-6 

1.08e-6 

1.27e-6 

1.27e-6 
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3.6.2  Astronomical  Requirements 

Time  simulations  are  performed  for  300  random  dynamics  within  the  uncertainty  range  (MonteCarlo 
analysis)  in  the  original  ESA  benchmark  simulator  (Bl)  described  in  Section  3.1  and  in  a  complementary 
benchmark  (B2),  both  developed  under  Matlab/Simulink.  The  latter  just  adds  to  Bl  a  low  frequency 
disturbance  force  at  plant  input  along  the  X-axis  (magnitude:  2  pN,  frequency:  oo  =  0.05  rad/sec)  in  order 
to  consider  disturbance  rejection  and  coupling  attenuation  at  low  frequencies.  In  each  simulation,  the 
criteria  appearing  in  Table  IV  are  computed  over  the  entire  simulation  time  (i.e.  5000  sec). 

In  order  to  characterize  the  minimum  performance  obtained,  the  worst  results  reached  by  every  controller 
are  presented  in  Table  V.  In  other  words,  for  each  controller,  the  greatest  value  over  the  300  uncertain 
cases  is  shown  for  the  maximum  and  the  standard  deviation  of  position  and  attitude  errors,  as  well  as  for 
maximum  actuator  commands,  in  all  axes.  Then,  it  is  possible  to  verify  whether  the  worst  performance 
still  respects  the  requirement.  The  bold  number  in  every  row  of  Table  V  represents  the  best  result  (best 
control  strategy)  for  every  particular  specification. 

Position  errors  (1,2, 3, 7 ,8,9  -Tabie  V) 

By  inspecting  Table  V,  it  is  found  that  the  performance  obtained  in  time  simulations  is  very  good 
concerning  position  accuracy,  since  the  requirements  are  easily  fulfilled  (an  improvement  of  two  orders  of 
magnitude  with  respect  to  the  specification  is  achieved  in  benchmark  Bl,  and  at  least  one  order  of 
magnitude  in  benchmark  B2)  for  maximum  and  standard  deviation  values.  The  non-diagonal  MIMO  QFT 
design  either  equals  or  slightly  improves  the  diagonal  MIMO  QFT.  Both  QFT  controllers  improve  the  H- 
inflnity  results  for  the  two  benchmarks. 

Attitude  errors  (4,5,6,10,11,12 -Tabie  V) 

The  specification  for  the  highest  attitude  error  is  harder  to  meet  mainly  because  of  the  effect  of  the  flexible 
modes.  Some  of  the  maximum  attitude  values  of  the  H-infinity  even  exceed  the  25.5  mas  required:  see 
benchmark  Bl  (4  -Table  V)  and  benchmark  B2  (4,5  -Table  V).  Again,  the  MIMO  QFT  methodologies 
improve  the  results  of  the  H-infinity  controller  in  the  six  attitude  error  cases  (4,5,6,10,1 1,12  -Table  V). 

Once  more,  the  non-diagonal  MIMO  QFT  either  equals  or  improves  the  diagonal  QFT  controller  results. 
The  greatest  differences  between  both  controllers  can  be  observed  at  the  Attitude  Error  along  the  Y-axis 
(5,11  -Table  V),  especially  for  benchmark  B2.  There,  the  non-diagonal  design  decreases  the  standard 
deviation  attitude  error  by  0.85  %  (1 1  -Table  V)  and  the  maximum  attitude  error  by  5.05  %  (5  -Table  V) 
with  respect  to  the  values  reached  by  the  diagonal  compensator.  These  improvements  could  turn  out  to  be 
relevant  to  the  astronomical  mission.  Their  achievement  is  due  to  the  fact  that  the  off-diagonal 
compensators  have  been  designed  to  minimize  the  coupling  at  low  frequencies,  which  are  principally 
stressed  in  the  second  benchmark. 


3.6.3  Engineering  Requirements 

Saturation  limits.  Actuator  commands  (13,14,15,16,17,18  -Table  V) 

As  can  be  seen  in  Table  V,  actuation  is  very  small  and  far  below  the  saturation  limits.  The  results  for  the 
three  controllers  remain  at  similar  values  (13,14,15,16,17,18  -Table  V). 

Open-loop  Bandwidth  Comparison 

The  open-loop  cross-over  frequency  results  of  the  six  SISO  loop  subsystems  are  shown  in  Table  VI.  These 
measures  correspond  to  the  smallest  frequencies  in  Hz  where  the  transfer  functions  of  the  open-loop  of 
each  SISO  subsystem  /?ii(s)  gii(s)  [without  the  coupling  elements  /^(s),  i  ^  j]  are  equal  to  0  dB.  The 
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minimum  performance  for  each  of  the  three  designs  has  been  established  as  the  minimum  bandwidth  value 
over  the  300  random  satellite  dynamics.  Obviously,  the  bandwidth  results  for  the  two  MIMO  QFT  designs 
coincide  since  their  diagonal  compensators  are  the  same.  A  value  of  0.01  Hz  is  considered  a  good 
compromise  choice  for  bandwidth.  Since  the  frequencies  of  the  first  flexible  modes  are  within  the  range 
[0.05,  0.5]  Hz,  the  open-loop  cross-over  frequencies  for  attitude  and  position  are  tuned  to  be  as  high  as 
possible  while  simultaneously  preventing  the  flexible  modes  from  disturbing  the  system  output 
performance. 


Table  VI.  Frequency  performance  with  the  three  controllers 


Requirement 

Non-diagonal 
MIMO  QFT 
Controller 

Diagonal  MIMO 
QFT  Controller 

H-infinity 

Controller 

Position  Bandwidth  X  (Hz) 

0.01  Hz 

0.0464 

0.0464 

0.0229 

Position  Bandwidth  Y  (Hz) 

0.01  Hz 

0.0472 

0.0472 

0.0206 

Position  Bandwidth  Z  (Hz) 

0.01  Hz 

0.0212 

0.0212 

0.0205 

Attitude  Bandwidth  X  (Hz) 

0.01  Hz 

0.00949 

0.00949 

0.0102 

Attitude  Bandwidth  Y  (Hz) 

0.01  Hz 

0.0102 

0.0102 

0.0102 

Attitude  Bandwidth  Z  (Hz) 

0.01  Hz 

0.00883 

0.00883 

0.0102 

max  r(/ffl|  (dB) 

CO 

<6  dB 

12.73  (Max) 

5.64  (Mean) 

11.51  (Max) 

5.55  (Mean) 

4.70  (Max) 

4.69  (Mean) 

max\s(jco^  (dB) 

CO 

<6  dB 

13.11  (Max) 

6.69  (Mean) 

12.05  (Max) 

6.60  (Mean) 

6.24  (Max) 

5.48  (Mean) 

For  the  position  transfer  functions,  the  three  controllers  exceed  the  0.01  Hz  recommendation  (two  and 
even  four  times  depending  on  the  controller  and  the  axis).  However,  the  flexible  modes  mostly  affect  the 
attitude  transfer  functions  and  do  not  impose  such  strong  constraints  on  the  position  transfer  functions. 
Consequently,  it  is  possible  to  go  over  0.01  Hz  for  the  position  loops,  as  is  proved  by  the  satisfying  time 
domain  results  in  Table  V. 


For  the  attitude  transfer  functions,  the  open-loop  cross-over  frequencies  are  around  0.01  Hz  for  the  H- 
inflnity  and  for  both  the  non-diagonal  and  diagonal  MIMO  QFT  designs. 


3.6.4  Control  Requirements 


Stability  Objectives 

Stability  and  performance  specifications  are  essentially  described  as  mathematical  expressions  ready  to  be 
used  during  the  design  process  of  the  controller.  These  expressions  usually  differ  from  one  control 
methodology  to  another  provided  they  are  based  on  distinct  approaches,  which  is  the  case  of  H-inflnity 
and  QFT-based  methodologies.  In  this  paper  the  stability  specifications  have  been  defined  in  two  different 
ways: 

a)  Stability  conditions  of  Section  2.6  for  MIMO  QFT  (Nyquist  criterion  for  sequential  methods). 

b)  Margins  on  max\T(jco\  and  max\s(jco\  for  H-infinity  (classical  criterion  for  MIMO  systems). 

(O  CO 


The  non-diagonal  and  the  diagonal  MIMO  QFT  controllers  fulfill  the  stability  conditions  for  sequential 
procedures  defined  in  Section  2.6.  The  H-infinity  compensator  fulfills  the  margins  of  T(s)  and  S(s )  defined 
in  Table  VI.  With  respect  to  this  classical  interpretation  of  robust  stability,  the  QFT  approaches  respect 
them  in  most  of  the  cases  (mean),  but  not  in  several  cases  (max).  This  is  due  to  the  fact  that  these 
interpretations  of  the  stability  margins  (which  are  indeed  a  margin  of  a  margin)  are  not  integrated  as  a 
design  specification  in  the  core  of  QFT  techniques. 
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Since  stability  and  performance  specifications  are  only  interpretations  of  the  functional  requirements 
(astronomical  and  engineering  requirements  -Table  IV-),  the  designer  should  be  aware  of  which  tradeoffs 
need  to  be  made.  Essentially,  the  interpretation  of  reality  in  terms  of  a  particular  theory  can  never  replace 
the  real  world  itself.  In  the  absence  of  the  real  system  implementation  or  a  suitable  prototype  to  be  used 
instead,  the  designer  must  manage  time  domain  simulations  in  order  to  verify  the  control  system  behavior 
[1].  That  interpretation  was  done  and  successfully  validated  in  the  previous  sections  of  the  paper. 

Additionally,  the  classical  margins  on  max\T(jco) I  and  max\s(jco\  are  stability  MIMO  margins,  but  they  do 

CO  CO 

not  include  phase  information.  This  fact  makes  them  sufficient,  but  not  necessary  conditions  and  could 
yield  very  conservative  controllers  in  some  situations.  Although  the  three  methods  are  robust  stable 
according  to  their  own  requirement  and  to  time  domain  simulations,  future  research  work  to  re-interpret 
both  types  of  robust  stability  conditions  and  margins  constitute  one  of  the  next  research  objectives. 

Controller  Complexity  and  Order 

The  number  of  operations  that  have  to  be  performed  per  sampling  period  may  place  restrictions  on  the 
compensator  design.  The  implementation  of  a  controller  based  on  the  state  space  representation  differs 
from  that  based  on  transfer  functions.  The  former  appears  to  have  a  common  denominator  for  every 
element  of  the  compensator  when  it  is  transformed  into  transfer  function  description  and  the  latter  does  not 
actually  need  it.  Indeed,  the  expression  of  each  control  signal  in  a  transfer  function  matrix  depends  on  its 
corresponding  row  of  the  compensator  matrix.  But  even  there,  common  denominators  are  not  needed.  The 
control  signal  Wi(v)  is  computed  as  the  sum  of  signals  generated  by  every  compensator  in  the  i- th  row. 

In  order  to  make  a  realistic  comparison  of  the  computational  cost  of  the  different  controllers  (non-diagonal 
MIMO  QFT,  H-infinity  and  diagonal  MIMO  QFT),  the  number  of  sums  and  multiplications  computed  in 
each  sample  at  the  final  implementation  are  analyzed.  Following  the  same  discretization  process,  the 
values  in  Table  VII  are  obtained.  The  compensator  matrix  of  the  H-infinity  design  expressed  in  transfer 
function  description  presents  36  elements  having  42nd  order.  The  diagonal  MIMO  QFT  design  consists  of 
six  diagonal  compensators  going  from  5th  to  14th  order.  The  non-diagonal  MIMO  QFT  design  consists  of 
eight  compensators  going  from  3rd  to  14th  order. 


Table  VII.  Computational  cost  per  sampling  for  the  three  controllers 


Controller 

Number  of  Multiplications 

Number  of  Sums 

Non-diagonal  MIMO  QFT 

130 

124 

H-infinity 

2994 

2988 

Diagonal  MIMO  QFT 

116 

110 

4.0  COMBINING  SWITCHING  &  ROBUST  QFT  CONTROL  STRATEGIES  TO 
IMPROVE  CLASSICAL  CONTROL[8] 

This  section  introduces  a  methodology  to  design  a  family  of  robust  controllers  able  to  go  beyond  the 
classical  linear  limitations.  Combining  robust  designs  and  stable  switching,  the  new  controllers  optimize 
the  time  response  of  the  system  by  fast  adaptation  of  the  controller  parameters  during  the  transient 
response  according  to  certain  rules  based  on  the  amplitude  of  the  error.  The  methodology  is  based  on  both 
a  new  graphical  stability  criterion  for  switching  linear  systems  and  the  robust  quantitative  feedback  theory 
(QFT). 
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Switching  control  has  demonstrated  to  be  an  efficient  tool  in  achieving  tight  performance  specifications  in 
control  systems  [15-16].  The  way  to  reach  this  enhancement  is  by  designing  various  parallel  controllers 
with  different  characteristics,  and  continuously  selecting  among  them  the  one  that  governs  the  system 
(Fig.  21).  Thus,  performance  specifications  that  are  not  achievable  by  a  simple  linear  controller,  as  the 
limitation  theory  predicts  [17-18],  can  be  attained  through  suitable  switching  rules. 

One  of  the  main  issues  in  switching  control  techniques  is  that  the  system  stability  is  not  assured  a  priori , 
even  if  the  switching  is  made  between  stable  controllers.  This  is  the  reason  why  most  of  the  current 
literature  about  switching  systems  is  still  devoted  to  stability  issues.  See  [19-20]  for  general  results  about 
stability  criteria  applied  in  some  particular  cases. 

4.2  Switching  &  robust  QFT  control 
4.2.1  Switching  Systems  Stability 


Fig.  21  shows  a  general  scheme  of  a  switching  system.  A  set  of  controllers  is  designed  and  a  supervisor 
selects  the  most  suitable  one,  depending  on  the  system  and  environment  parameters. 


Controllers 


Fig.  21  Switching  control  scheme. 

One  of  the  main  difficulties  found  when  switching  techniques  are  applied  is  that,  in  general,  the  system 
stability  is  not  assured,  even  if  switching  is  made  between  stable  controllers.  Some  extra  conditions  must 
be  met.  In  particular,  it  has  been  proved  that  a  system 

x(t)  =  A A(Q  e  A  =  {Alv..,Aj,  A.  Hurwitz,  (71) 


with  arbitrary  switching  within  the  set  of  matrices  A  is  exponentially  stable  if  and  only  if  there  exists  a 
common  Lyapunov  function  (CLF)  for  all  A*  in  the  set  A  [22].  It  has  also  been  proved  that  the  existence  of 
a  common  quadratic  Lyapunov  function  (CQLF)  is  a  sufficient  condition  for  exponential  stability  [23].  In 
this  context,  the  main  issue  in  linear  switching  systems  is  finding  conditions  under  which  the  existence  of 
a  CQLF  is  assured.  In  particular,  it  has  been  proved  that  the  circle  criterion  provides  necessary  and 
sufficient  conditions  for  the  existence  of  a  CQLF  for  two  systems  in  companion  form  [24-26],  that  is,  the 
systems 


(72) 


(73) 
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both  Hurwitz,  with 
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have  a  CQLF  if  and  only  if 

1  +  Re{Ar  (si  -  A)-1  g}  >  0,  s  =  yco,  for  all  frequency  cq.  (75) 


This  paper  considers  stability  for  arbitrary  switching  between  two  closed-loop  systems  with  transfer 
functions  Ti(s)  =  Zi(s)/[1+Zi(s)]  and  T2(s)  =  Z2(s)/[1+Z2(s)],  both  stable,  where  Zi(s)  =  P(s)Gi(s)  and  L2{s) 
=  P(s)G2(s)  are  proper,  have  the  same  number  of  poles,  and  the  same  number  of  zeros.  In  this  case,  the 
effect  of  switching  is  to  change  the  gain  and  the  position  of  poles  and  zeros.  The  open  loop  transfer 
functions  of  both  systems  are 


Lx(s)  = 


b^+...  +  b0 

S  +  &n_\S  +...  +  £Zq 


N(s) 

D(s) 


and 


(76) 


(&„_!  +A b^s”  1  +...  +  (60  +A b0)  _  N(s )  +  AN (5) 

sn  +  +  A ani ) s "  1  + ...  +  (a0  +  Aa0)  D(s)  +  A D(s^ 


Note :  For  the  sake  of  clarity  in  the  subsequent  analysis,  a  general  expression  has  been  used  where  the 
order  or  the  numerator  is  one  less  than  that  of  the  denominator.  If  it  were  not  the  case,  then  the  coefficients 
bn.  1,  etc,  would  be  zero. 

The  closed-loop  transfer  functions  are 


T](s)  = 


N(s) 


LM  _ 

1  +  Z,  (x)  D(s)  +  N(s) 


and 


T2(s)  = 


L2  (s) 


N(s)  +  AN(s) 


l  +  L2(s)  D(s)  +  AD(s)  +  N(s)  +  AN(s ) 
where  the  characteristic  equations  are 

—  sn  +  &n_iSn  1  + ...  +  c^s  +  and 


D{s)  +  AD(s)  +  N(s)  +  AN  (v)  —  sn  +  (en_x  +  A  en_x)sn  + ...  +  (e0  +  Ae0), 


with  e,  =  a,  +b;  and  Ae,  =  A  a,  +  A  b,. 


1  1  1 


1  1 


(78) 


(79) 


(80) 

(81) 
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Using  these  expressions  for  the  coefficients  eL  and  Aeh  the  matrices  A,  g,  and  A  are  defined.  Now  the  circle 
criterion  is  applied  to  guarantee  stability  under  arbitrary  switching.  As 

01-  A)-1  =  sn  +  en_xsn~x  +  ...  +  e1s  +  e0  and  Arg  =  A en_xsn~l  +  ...+  Aexs  +  Ae0,  (75)  becomes 


1  +  Re< 


A.cn_xsn  + . . .  +  Acxs  +  A^q 
sn  +  cn_xsn  + ...  +  cxs  +  Cq 


>  0 ,  S  —  j CO,  for  all  frequency 


co. 


(82) 


After  some  simple  manipulation, 


Res  1  + 


A £n xs  + ...  +  A  cxs  +  Ae0 

S  +  &n_\S  +  ...  +  CXS  + 


=  Re 


f  N(s)  +  A N(s)  +  D(s)  +  ADO) 


r 


N(s)  +  D(s ) 


=  Re 


N(s)  +  ANjs )  +  £>(5)  +  ADCs) 
N(s)  +  AN(s) 


N(s)  +  D(s) 


rN(s)  +  AN(s)^ 
N(s) 


=  Re 


1  +  4(0 


1  +  4(0 


D(s)  +  A£>4) 

Z)(0 


(83) 


the  condition  can  be  expressed  in  the  following  form: 

>0,  5  =  j  co,  for  all  frequency  co. 


Res 


fi+44) 

(  z>(s)+ad(s)Y[ 

[i+400 

l  m  Jj 

(84) 


As  this  formulation  of  the  circle  criterion  is  applied  to  open-loop  transfer  functions,  it  is  enough  to  check  it 
at  positive  frequencies  due  to  symmetry.  Condition  (84)  is  then  equivalent  to 


arg  {1  +  4  C/g>)}  -  arg  {1  +  Lx  (yco)}  +  arg- 


D(j<J>)  +  APC/eo) 
D(M 


<  —  for  all  (0  >  0. 

2 


(85) 


Let  us  denote 


cpi2(©)[deg]  =  |arg{l  +  4(yco)}  -  arg{l  +  4(yco)}|, 


(86) 


a(co)  [deg]  = 


D(M  +  ADC/co)  1 

D(M  J 


(87) 


Using  the  triangle  inequality,  a  sufficient  condition  for  (85)  is 

cp12(co)  <  90  —  a(co)  deg  for  all  co  >  0.  (88) 

Then,  the  criterion  can  be  applied  graphically  in  both  the  Nyquist  and  the  Nichols  diagrams.  In  the  first 
case,  the  criterion  may  be  expressed  by  saying  that  £i(/oo)  and  L2(j(n)  must  be  inside  of  an  arc  of 
[90  -  a(oo)]  deg  around  the  point  (-1,0)  at  each  frequency.  In  the  Nichols  diagram  a  condition  over  angles 
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is  more  easily  checked:  plotting  the  frequency  response  of  [1  +  ii(/’co)]  and  [1  +  £2(/cq)],  the  distance 
912(00)  on  the  horizontal  axis  at  each  frequency  must  be  less  than  [90  -  a(oo)]  deg. 

It  must  be  noted  that  the  function  a(oo)  can  be  expressed  in  the  following  form: 


a(co)  = 


arg< 


DQco)  +  ADQ'co)] 

.  J 


arg- 


YlUa  +  Pj+APj) 

7=1 


Ylu^  +  Pj) 


j= 1 


(89) 


Z  ar§  t/'©  +  Pj  +  APj  1  -  ar§  (7©  +  Pj } 


7=1 


so  it  may  be  considered  as  a  measurement  of  the  controller  poles  change,  as  Fig.  22  shows.  Consequently, 
the  larger  the  movement  made  by  the  poles,  the  bigger  the  conservativeness  introduced  by  the  triangle 
inequality  at  (88).  Note  that  for  each  frequency  there  is  a  different  angle  a(oy). 


Im 


- - -  ^  -V?'' 

-p2  a2(coz)\  ^  ^  C  I 

„  "  / x  /  /  / 

«,+4 


X 


M 


,  -P'  -(P.+Api) 

,  ,'u3(md 

-(P,+Ap,)  y 


(0,0)  Re 


X 

-p. 


Fig.  22  a (s)  for  a  system  with  three  switching  poles.  a(co*)  =  a^co,)  +  a2(co/)  +  astco,). 


At  this  point  two  questions  arise.  Firstly,  the  criterion  presented  above  can  be  applied  to  switching 
between  two  isolated  controllers  with  the  same  structure.  However,  it  is  possible  that  the  designer  wants  to 
do  switching  among  more  than  one  system,  or  even  among  an  infinite  number  of  systems,  which  can  also 
be  considered  as  a  linear  parameter  varying  (LPV)  system,  where  the  controller  varies  continuously. 
Secondly,  real  systems  present  uncertainty,  so  the  criterion  must  be  modified  in  some  way  to  take  the 
uncertainty  into  account.  The  next  discussion  undertakes  both  issues. 

If  the  switching  is  made  among  a  set  of  controllers,  the  criterion  has  to  be  accomplished  by  every  pair  of 
them.  Checking  this  condition  may  be  an  impossible  task  if  there  is  more  than  one  pole  moving,  because 
the  angle  a  is  different  for  each  pair  of  controllers.  For  this  reason,  if  we  are  interested  in  a  controller 
whose  parameters  change  continuously  with  the  error,  we  will  permit  only  variable  gain  and  zeros.  Then, 
the  angle  a(co)  is  null  for  every  frequency,  and  the  only  condition  to  satisfy  is  that  the  angle  between  any 
two  possible  systems  Li(j co)  and  the  critical  point  (-1,0)  is  less  than  90  deg.  Moreover,  under  this  premise 
the  conservativeness  introduced  in  (88)  vanishes.  The  condition  can  be  checked  graphically  with  a  grid  of 
the  possible  open-loop  systems  that  the  controller  variation  can  generate,  as  Fig.  23a  shows.  The 
maximum  angle  cpi2(co)  must  be  contained  in  a  90  deg  arc  from  (-1,0).  In  the  Nichols  Plot,  the  way  to 
apply  the  criterion  is  to  draw  the  grid  of  possible  1  +£z(/a>)  systems,  and  check  that  the  maximum 
horizontal  distance  is  less  than  90  deg.  Fig.  23b  illustrate  the  criterion  in  the  Nichols  Plot. 
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Using  similar  arguments,  it  is  also  easy  to  deal  with  uncertainty.  For  an  uncertain  system,  the  template 
3P(/'co)  is  the  area  of  the  possible  plants  within  the  uncertainty  at  the  frequency  oo.  If  this  system  is 
governed  by  a  switching  controller,  each  point  of  the  template  can  change  its  position  due  to  switching. 
From  this  point  of  view,  switching  can  be  considered  as  a  mechanism  that  modifies  the  position  and  shape 
of  the  templates  of  [1  +  Zz(/co)].  To  be  sure  that  the  switching  is  stable,  the  above  criterion  must  be  applied 
to  the  whole  template. 


dB 


90  deg 


Fig.  23  Criterion  for  continuous  switching:  a)  Complex  plane,  b)  Nichols  plot. 

It  has  been  traditionally  considered  in  control  theory  that  uncertainty  changes  the  plant  slowly  in 
comparison  with  the  system  dynamics.  If  the  switching  laws  depend  on  the  state  of  the  system,  then  the 
switching  is  much  faster  than  changes  due  to  uncertainty.  Consequently,  for  each  point  of  the  departure 
template  there  is  only  one  corresponding  point  in  the  arrival  template.  Furthermore,  it  can  be  assumed  that 
uncertainty  does  not  affect  the  angle  a. 


Then  the  Nichols  Chart  is  a  very  clear  diagram  to  test  the  stability  of  the  uncertain  switching  system.  Fig. 
24  shows  the  templates  of  [1 +/,*•(/©)]  and  the  application  of  the  method.  If  during  the  displacement  of  each 
point  of  the  first  template  to  its  corresponding  point  of  the  second  one,  the  maximum  horizontal  distance 
between  any  two  points  of  this  path  is  less  that  90  deg,  the  stability  condition  is  satisfied  at  that  particular 
frequency.  Although  it  is  a  laborious  task  to  check  each  point  of  the  template  at  each  frequency,  usually  it 
is  not  necessary  because  the  whole  set  of  templates  are  much  closer  together  than  the  critical  distance. 


Fig.  24  Stability  criterion  on  the  Nichols  Chart 
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4.2.2  Methodology 

Quantitative  feedback  theory  [1]  has  demonstrated  to  be  an  excellent  tool  dealing  with  the  compromises 
between  several,  often  conflicting,  control  specifications.  Its  transparent  design  process  allows  the 
designer  to  consider  all  of  them  simultaneously  and  in  the  same  plot.  The  QFT  philosophy  permits  to 
design  a  controller  that  satisfies  the  set  of  performance  specifications  with  every  plant  within  the 
uncertainty,  using  the  minimum  amount  of  feedback. 

However,  as  any  linear  methodology,  QFT  produces  linear  controllers,  and  consequently  it  is  subject  to 
the  performance  limitations  of  linear  systems.  Briefly  speaking,  the  way  in  which  a  system  can  improve  its 
performance  is  limited  by  the  reduction  of  the  stability  margins.  Although  those  limitations  are  usually 
expressed  in  the  frequency  domain,  they  have  consequences  in  the  time  response  of  the  system. 

As  limitations  are  due  to  the  linear  character  of  the  system,  it  seems  that  the  key  is  to  use  nonlinear 
controllers.  Some  work  has  been  done  in  that  respect,  and  some  elegant  solutions  have  been  performed 
[27-28].  However,  the  use  of  nonlinear  elements  may  imply  some  problems,  like  instability  or  limit  cycles, 
as  well  as  the  lack  of  simple  design  tools. 

Another  solution  to  overcome  the  linear  limitations  is  to  prioritize  some  specifications  over  others 
according  to  the  state  of  the  system  at  each  time.  In  other  words,  do  switching  to  select  the  appropriate 
controller,  keeping  always  the  linear  character.  A  particular  way  to  do  this  is  to  use  the  error  amplitude  as 
the  switching  signal.  Then,  when  the  output  is  far  from  the  reference,  the  system  needs  to  be  more  stable, 
and  also  faster,  but  precision  is  not  so  necessary.  Conversely,  when  there  is  little  error,  some  amount  of 
stability  margin  can  be  sacrificed  in  order  to  increase  the  low  frequency  gain,  and  therefore  the  precision 
and  the  disturbance  rejection.  These  ideas,  combined  with  the  QFT  method,  lead  to  a  new  design 
procedure,  listed  in  the  next  four  steps: 

Step  1 :  Obtain  a  preliminary  linear  controller  for  the  system,  usually  by  applying  QFT:  represent  the 
parametric  and/or  non-parametric  uncertainty  with  the  templates,  define  the  frequency  domain 
specifications,  generate  the  QFT-bounds  and  design  a  linear  controller  by  loop-shaping. 


Step  2:  The  preliminary  QFT  controller  is  the  starting  point  to  design  two  extreme  controllers  with  the 
same  structure,  where  gain  and  zeros  can  vary  freely,  but  poles  stand  still.  The  characteristics  of  these  two 
controllers  must  be  related  with  the  error  amplitude.  As  Boris  Lurie  explains  [16],  when  the  error  is  large 
the  bandwidth  must  increase  to  get  fast  response,  but  the  loop  gain  does  not  need  to  be  high.  And  for  small 
errors,  the  bandwidth  is  reduced  to  avoid  the  effects  of  noise,  while  the  low  frequency  gain  is  increased  to 
minimize  the  jitter  and  the  tracking  error. 

In  terms  of  loop-shaping,  these  rules  can  be  considered  as:  1)  for  small  errors,  increase  gain  and  move 
further  away  the  zeros,  and  2)  for  high  errors,  decrease  gain  and  bring  nearer  the  zeros.  Apart  from  this, 
reasonable  stability  margins  must  be  maintained,  although  they  could  be  considerably  reduced  for  the 
small  error  situation. 

Step  3:  The  robustness  of  the  extreme  designs  guarantees  that  both  linear  systems  are  stable  for  every  plant 
within  the  uncertainty.  However,  it  is  necessary  to  apply  the  criterion  presented  in  the  previous  section 
[see  equation  (88)  with  a  =  0]  to  assure  that  the  switching  between  both  controllers  is  also  stable. 

One  advantage  of  this  graphical  criterion  is  that  it  gives  information  about  the  frequencies  where 
conditions  are  not  satisfied,  so  the  designer  can  go  back  to  Step  2  and  change  the  extreme  controllers  in 
this  region. 
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Step  4\  Select  the  switching  function.  Simulations  of  the  system  governed  by  each  controller  are 
performed.  Looking  at  the  time  response  of  the  designs,  the  function  that  relates  the  error  amplitude  with 
the  position  of  the  controller  parameters  is  designed  (see  for  example  Fig.  27). 

4.3  Application:  Remotely  Controlled  Reconnaissance  Vehicle 

This  section  presents  a  simple  example  to  illustrate  the  new  methodology  introduced  in  section  4.2.  It 
consists  in  a  remotely  controlled  reconnaissance  vehicle  (see  [29]  for  details). 

The  plant  to  be  controlled  is: 


P(s)  = 


{s2  +  a{s  +  a0 


where  a{  e  [1.8, 2.2],  and  a0  e  [3.6, 4.4] 


Step  1 

The  controller  designed  by  Dorf  (see  [29])  is: 


(90) 


G0(s)  = 


10(5+2) 

“ML 


(91) 


Step  2 

Now,  starting  from  G0(s),  the  extreme  controllers  Gi(s)  and  G2(s ),  for  low  and  high  errors  respectively,  are 
designed  observing  the  guidelines  exposed  in  section  4.2.  The  elements  to  be  varied  are  the  gain  and  a 
zero.  The  new  controllers  are 


G{(s)  = 


15(5  +  2.3) 


and 


G2(5)  = 


8  (^  +  0.5) 

“mT 


(92) 


(93) 


and  their  frequency  domain  characteristics  can  be  viewed  in  the  Bode  plot  of  Fig.  25.  As  is  appreciated, 
Gi(s)  presents  a  higher  low-frequency  gain  and  a  dominant  pole,  while  in  G2(s )  the  low-frequency  gain  is 
lower  and  the  zero  is  the  dominant  one.  Their  robust  stability  margins  are  checked  with  QFT  tools. 


Fig.  25  Bode  plot  of  the  three  controllers 
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Step  3 

Fig.  26  shows  the  templates  T[  1  +  L\(jco)\  and  T[  1  +  L2(jco )]  at  various  representative  frequencies  [1,  2,  3, 
5]  rad/sec,  and  the  Nichols  plot  of  each  nominal  function  [1  +  L;0(jco)\.  As  there  is  no  variation  in  any  pole 
of  the  controller,  the  angle  a  is  cero  for  all  frequencies,  and  the  only  thing  to  check  is  that,  in  the  path  from 
each  point  of  the  first  template  to  its  corresponding  point  of  the  second  template,  the  maximum  horizontal 
distance  between  two  points  is  not  higher  than  90  deg  [see  (88)  with  a  =  0]. 


Fig.  26  Stability  study  of  the  switched  systems  on  the  Nichols  chart 


Step  4 

The  switching  controller  Gswi{s )  presents  the  following  expression: 

(l5-7q(s  +  2.3-1.8fc) 

swA)~  (^  +  1) 


(94) 


where  the  parameter  k  is  given  by  a  function  Z  — >  [0,  1]  of  the  error  signal.  In  order  to  reduce  possible 

impulse  effects,  a  smooth  function  (95)  has  been  selected  instead  of  a  relay-type  or  saturation-type 
function, 


f 

k  =  1  -  exp 

v 


*(02> 

0.15, 


(95) 


where  the  values  of  the  parameters  have  been  adjusted  for  optimal  performance  by  simulation.  Some 
further  research  could  be  devoted  to  finding  the  optimum  shape  for  the  switching  functions.  The  shape  of 
the  switching  function  is  shown  in  Fig.  27. 
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Fig.  27.  Switching  function 


Validation 

The  compensators  are  verified  in  the  time  domain  for  the  most  representative  plants  selected  from  the  set 
of  uncertain  plants.  Fig.  28  shows  the  time  response  of  the  control  system  (nominal  plant  [29])  to  a  step 
input  reference  tracking  (t  =  1  sec)  and  a  step  disturbance  (t  =  15  sec).  The  response  of  the  controller 
designed  with  the  new  methodology  (switching  controller,  Gswi)  combine  the  best  characteristics  of  the 
extreme  controllers  ( G\  and  G2),  improving  the  response  of  the  original  fixed  controller  (G0). 


Fig.  28  Response  to  a  step  input  reference  tracking  and  a  step  disturbance 


5.0  CONCLUSIONS 

Since  the  very  first  ideas  suggested  by  Horowitz  in  1959  until  now,  the  Quantitative  Feedback  Theory 
(QFT)  has  been  successfully  applied  to  many  control  systems:  linear  and  non-linear,  stable  and  unstable, 
SISO  and  MIMO,  minimum  and  non-minimum  phase,  with  time-delay,  with  lumped  and  distributed 
parameters,  multi-loop,  etc  [1].  The  method  searches  for  the  controller  that  guarantees  the  achievement  of 
the  required  performance  specifications  for  every  plant  within  the  existing  model  uncertainty.  QFT 
highlights  the  trade-off  (< quantification )  among  the  simplicity  of  the  controller  structure,  the  minimization 
of  the  ‘cost  of  feedback’,  the  quantified  model  uncertainty  and  the  achievement  of  the  desired  performance 
specifications  at  every  frequency  of  interest. 

The  first  part  of  the  paper  summarized  a  methodology  to  design  sequential  non-diagonal  QFT  controllers 
for  multi-input-multi-output  MIMO  systems  with  uncertainty.  The  second  part  demonstrated  the  feasibility 
of  that  methodology  to  control  UAVs,  in  particular  the  position  and  attitude  control  of  a  6x6  MIMO 
Darwin-type  spacecraft  with  large  flexible  appendages.  The  third  part  of  the  paper  introduced  a  new 
practical  methodology  to  design  robust  controllers  that  work  under  a  switching  mechanism.  The  method  is 
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capable  of  optimising  performance  and  stability  simultaneously,  going  beyond  the  classical  linear 
limitations  and  giving  a  solution  for  the  well-known  robustness-performance  trade-off.  Based  on  the 
frequency  domain  approach,  the  method  combines  a  graphical  stability  criterion  for  switching  linear 
systems  and  the  robust  quantitative  feedback  theory. 
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